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INTRODUCTION 



Les réseaux hydrographiques constituent les principaux vecteurs d'évacuation de la 
matière à la surface des continents. Sur un plan géomorphologique, ils sont le siège d'une 
érosion active due aux écoulements d'eau concentrés qui permettent un transport rapide 
de sédiments sur de longues distances. Us s'opposent ainsi aux versants, où l'érosion 
diffuse domine et engendrent une érosion différentielle à l'origine de la formation des 
talwegs et des vallées. Ce phénomène conduit à long terme à un découpage de la 
topographie suivant une grande variété de longueurs d'onde. Ils constituent ainsi un des 
principaux agents de structuration du relief. Sur un plan hydrologique, les réseaux 
hydrographiques en drainant les écoulements de surface et de subsurface contrôlent 
l'humidité générale d'un bassin versant et accélèrent les transferts de volumes d'eau au 
sein des bassins. Ils constituent à ce titre l'un des paramètres importants dans la 
génération des pics de crue et plus généralement l'une des composantes principales de la 
réponse hydrologiques d'un bassin versant aux précipitations. 

L'efficacité des réseaux hydrographiques à évacuer de la matière (eau ou sédiments) 
dépend étroitement d'un part de l'énergie cinétique des écoulements dont ils sont le siège 
et d'autre part de leur répartition spatiale au sein de la surface topographique. La densité 
des réseaux, en imposant une longueur moyenne d'interfluve, fixe en effet une limite 
supérieure pour les temps de transfert "lents" et par là même des taux d'évacuation. La 
dynamique de l'évolution des reliefs est donc conditionnée par l'organisation spatiale des 
réseaux hydrographiques. 

D'autre part, les réseaux hydrographiques traduisent la réponse des processus 
d'érosion à une augmentation de l'énergie potentielle de la surface topographique induite 
par des mouvements tectoniques ou eustatiques. Ces processus d'érosion répondent ainsi 
à des critères d'équilibre et tendent à faire évoluer la topographie vers une configuration 
stable fixée par des conditions aux limites données. La géométrie des systèmes témoigne 
dans certains cas de la nature des processus physiques qui les engendrent. Nous pouvons 
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donc espérer, à travers l'organisation géométrique des réseaux hydrographiques, 
retrouver les lois physiques d'évolution et d'organisation de la topographie. 

La compréhension de l'organisation des réseaux hydrographiques est donc un enjeu 
important tant sur le plan fondamental dans l'étude de la dynamique de l'évolution de la 
surface des continents, que sur le plan pratique dans une perspective d'application aux 
prévisions hydrologiques. 

L'étude de la géométrie des réseaux hydrographique vise à établir: (i) les règles de 
connexion des arborescenceîhydrographiques, (ii) les distributionîde longueur de cours 
d'eau ou de paramètres topographiques tels que l'aire drainée, (iii) les corrélations entre 
les différents paramètres géométriques. Ces relations doivent autant que possible être 
indépendantes de l'échelle spatiale pour définir des propriétés générales des systèmes 
hydrographiques. Un certain nombre de lois géométriques ont pu ainsi être observées et 
utilisées comme critère de validation dans des modèles de construction de réseaux 
hydrographiques (e.g. Horton, 1945; Hack, 1957). Cependant ces lois apparaissent 
insuffisantes pour construire un modèle physique d'organisation. En effet, elles ne 
caractérisent pas les propriétés intrinsèques de ces réseaux mais plutôt celle des 
arborescences dans un sens très générale. Il est donc nécessaire d'affiner l'analyse 
géométrique et statistique de manière à en extraire les caractéristiques propres aux réseaux 
hydrographiques. C'est l'objet du travail présenté dans le premier chapitre de ce mémoire 
où une étude statistique détaillée de deux réseaux hydrographiques bretons (la Vilaine et 
le Blavet) a été réalisée. 

L'analyse des réseaux hydrographiques actuels se réfère cependant à un état figé du 
système et ne permet pas de remonter à sa dynamique d'élaboration au cours du temps. 
La caractérisation géométrique des réseaux est donc insuffisante pour une description 
complète de leur organisation. Pour appréhender la dynamique du système, il est donc 
nécessaire de construire des modèles de croissance des réseaux hydrographiques. Une 
des bases de ces modèles est de considérer que l'extension des réseaux hydrographiques 
est conditionnée par la topographie, cependant la topographie est elle même transformée 
par l'action érosive des cours d'eau. Il est donc important de tenir compte du lien 
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étroit entre les formes du relief et le réseau hydrographique via les processus d'érosion. 
Les modèles proposés dans la littérature (Wilgoose et al., 1991abc; Howard, 1994) sont 
construit sur une approche détaillée des processus et mettent en oeuvre un nombre 
important de paramètres. Nous avons préféré une approche basée sur une modélisation 
simple, dans laquelle nous distinguons deux types de phénomènes: (i) les mécanismes 
d'érosion dans le système hydrographique, qui sont générateurs d'une redistribution des 
volumes sur de grandes distances; (ii) les mécanismes d'érosion hors des systèmes 
hydrographiques (sur les "versants"), qui sont à l'origine d'une érosion et de 
déplacements locaux et diffus. Nous nous hmitons à quatre paramètres de premier ordre, 
que sont le débit critique de formation des cours d'eau, la non linéarité des processus 
d'érosion sur le versant et dans le réseau, le coefficient de diffusion sur le versant et la 
capacité de transport des cours d'eau. Dans le deuxième chapitre de ce mémoire, nous 
présentons ce modèle et son utilisation pour l'étude de l'influence des paramètres sur 
l'évolution du relief et des réseaux hydrographiques. 

Les deux premières parties de ce travail mettent en évidence que la définition des 
limites amonts des cours d'eau constitue un paramètre déterminant dans l'organisation 
des structures hydrographiques. Elle détermine en effet la position géographique de la 
transition versant-réseau. Dans les régions humides telles que la Bretagne, les sources et 
les abords du réseau hydrographique se présentent sous la forme de zones saturées en 
eau. Il est donc difficile de localiser précisément cette transition. L'objet de la troisième 
partie de ce travail est de proposer une approche expérimentale visant à définir les 
paramètres topographiques permettant une localisation de la transition versant-réseau dans 
le cas de régions humides. Nous avons ainsi cherché à déterminer dans quelle mesure 
l'extension des zones saturées reflète celle du réseau hydrographique. 
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CHAP. 1: ETUDE STATISTIQUE DE L'ORGANISATION DES 
RÉSEAUX HYDROGRAPHIQUES 



EXEMPLE DE LA VILAINE ET DU BLAVET 



1. Objectifs et limites de l'étude quantitative 

Contexte de V étude 

Le réseau hydrographique est à la fois le résultat et le principal acteur du remodelage 
de la topographie. De ce fait, il représente la signature d'un certain nombre de processus 
de transfert d'eau et d'érosion. Aussi est-il raisonnable d'envisager d'obtenir des 
informations sur la nature des phénomènes physiques liés à l'érosion et les limites 
spatiales de leur action à travers l'étude de la géométrie de ce réseau. 

Dès les années 50 les géomorphologues ont quantifié la géométrie des réseaux afin 
d'établir des lois statistiques propres aux systèmes hydrographiques (Horton, 1945; 
Morisawa, 1962,1964; Melton, 1958; Hack, 1957; Langbein, 1947). Si ces dernières 
correspondent à des caractéristiques plus ou moins générales de l'organisation 
géométrique des rivières dans un paysage, leur interprétation physique reste cependant 
très limitée. Certains travaux se sont également attachés à montrer et à définir la nature 
fractale (géométrie invariante d'échelle) du réseau hydrographique pris dans sa globalité 
et du parcours filaire d'une rivière (Hjelmfelt, 1988; Tarboton et al., 1988, 1990; La 
Barbera et Rosso, 1989, 1990; Rosso et al., 1991; Robert et Roy, 1990; Nikora, 1991, 
1994). Si la plupart des études adhèrent à la thèse d'un réseau fractal, il faut cependant 
rester prudent sur la représentativité des résultats obtenus jusqu'à présent (Nikora, 1991). 

Ces relations empiriques et fractales, à défaut de pouvoir donner des indices sur la 
physique qui les engendre, sont considérées comme des propriétés spécifiques aux 
réseaux hydrographiques. A ce titre elles sont utilisées comme critères de validation dans 
les modèles de simulation de construction de réseaux. Plusieurs modèles basés sur des 
hypothèses différentes physiques (Leopold et Langbein, 1962; Meakin et al., 1991; Roth 
et al., 1989; Howard, 1990; WiUgoose et al., 1991abc; Rodriguez Iturbe et al, 1992; 
Rinaldo et al., 1992; Rigon et al.; 1993) aboutissent dans leur grande majorité à des 
réseaux théoriques vérifiant les lois empiriques des réseaux naturels. Donc deux 
constatations s'imposent: (i) la relation entre la géométrie des réseaux et la physique qui 
l'engendre n'est pas univoque, (ii) il est toujours nécessaire de rechercher des lois 
statistiques ou géométriques, propres aux réseaux hydrographiques naturels, ou propres 
à des processus types (par exemple aléatoires ou déterministes). 
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Le travail présenté dans ce chapitre vise à obtenir des données sur la géométrie de 
réseaux hydrographiques d'une superficie d'un voir deux ordres de grandeur supérieurs 
(10 000 km2) aux travaux antérieurs. Nous pourrons ainsi, grâce à une meilleure 
représentativité des données, contribuer à mieux définir (i) les lois statistiques et 
l'organisation caractéristique des réseaux, (ii) les échelles pour lesquelles ces lois sont 
valables, (iii) les processus physiques sous-jacents. 

1.2. Choix des bassins versants de la Vilaine et du Blavet 

Pour des raisons d'intérêt régional et de disponibilité des informations, notre attention 
s'est portée sur deux réseaux hydrographiques bretons. La Bretagne ne présentant que 
très peu de bassins versants supérieurs à 2000 km^, la Vilaine et le Blavet furent 
sélectionnés en fonction de leur position dans la hiérarchie en terme d'aire drainée des 
grands bassins versants bretons. 

1.3. Intérêt du couplage réseau vectorisé-MNT 

Les données sont de deux types: 

- une topographie digitalisée (Modèle Numérique de Terrain: MNT) de la 
Bretagne au pas de 250 m qui englobe les bassins versants de la Vilaine et du 
Blavet, acquis grâce au GDR (Groupe De Recherche CNRS "Agriculture et 
Environnement"), 

- les réseaux hydrographiques de la Vilaine et du Blavet vectorisés sur la base 
de relevés de l'Institut Géographique National (IGN) au 1/100 000, fournis par 
l'agence de bassin Loire-Bretagne. 

Ces deux fichiers offrent le double intérêt de couvrir une superficie suffisamment 
grande, et d'avoir une résolution spatiale acceptable. Ils permettent ainsi d'aborder les 
problèmes d'échelle des structures géométriques des bassins versants, sur plusieurs 
ordres de grandeur. 

Les travaux précédents dans ce domaine font référence à des mesures de longueurs et 
de superficies calculées, soit sur des documents cartographiés (Hack, 1957; Melton, 
1958; Smart, 1969b, 1972; Abrahams, 1972; Jarvis, 1976; Robert et Roy, 1990), soit à 
partir de réseaux théoriques de drainage issus de MNT (Tarboton et al., 1988, 1991; 
Depratère et Moniod, 1991; Moussa, 1991; Andrle, 1992). L'utilisation simultanée d'un 
fichier vectorisé et d'un MNT permet d'allier les avantages des deux sources 
d'informations, à savoir: 

- une détermination indépendante du réseau hydrographique vrai (et non pas 
d'un pseudo réseau de drainage), 

- le traitement de très grandes surfaces, par rapport aux études précédentes. 
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1.4. Limites de l'étude 

La Bretagne présente un relief assez mou avec une longueur moyenne d'interfluve de 
l'ordre de 1 km. Le profil d'un versant n'est donc décrit que par 2 à 4 pixels de 250 m de 
côté. Il sera donc difficile, voire irréaliste, dans le cadre de ce travail, d'aborder l'étude de 
la géométrie des objets topographiques à des échelles inférieures ou égales à celles des 
versants. C'est pourquoi nous ne pourrons pas effectuer la statistique, ô combien 
intéressante, des paramètres topographiques susceptibles de contrôler l'existence d'un 
cours d'eau: tels que la convexité des talwegs (Depratère et Moniod, 1991), ou les 
paramètres dérivés de la pente locale (Dunne, 1980; Kirkby, 1988; Montgomery et 
Dietrich, 1988, 1992, 1994). Ce chapitre portera essentiellement sur deux paramètres 
géométriques: les longueurs des cours d'eau et l'aire drainée. Pour simplifier les calculs, 
l'aire drainée est calculée dans le plan d'échantillonnage du MNT (et non pas dans un 
espace à 3 dimensions). La Bretagne ne présentant pas de relief très prononcé, l'écart 
entre la valeur réelle de l'aire drainée et la surface calculée dans un plan, sera considéré 
comme négligeable devant celui induit par un relevé de l'altitude tous les 250 m. 

2. Bassins versants de la Vilaine et du Blavet 

2.1. Contexte géologique 

La description géologique qui suit est loin d'être exhaustive. Elle a pour simple but de 
décrire les grand traits géologiques du Massif Armoricain afin d'établir les principales 
caractéristiques des bassins versants de la Vilaine et du Blavet. Le lecteur trouvera une 
description plus détaillée dans le travail de synthèse de Le Corre et al. (1991) dont 
s'inspirent très largement les paragraphes qui suivent. 

2.1.1. Le Massif Armoricain 
Le Massif Armoricain est pour l'essentiel, le résultat de deux évolutions 
géodynamiques successives: 

-l'histoire précambrienne marquée par l'orogénèse cadomienne (620 à 540 Ma); 
-l'histoire paléozoïque marquée par l'orogénèse hercynienne (540 à 290 Ma). 

Il est limité à l'Est par les sédiments mésozoïques du bassin de Paris, et vers l'Ouest 
par les côtes de la Manche et de l'Atlantique (Fig. 1.1). Ces limites ne traduisent 
cependant pas la véritable étendue de ce massif qui se prolonge en réalité sous le Bassin 
Parisien et sur le plateau continental péri-armoricain. 

Il est composé de quatre grands domaines (Nord-Armoricain, Médio- Armoricain, 
Sud- Armoricain, Nord-Ouest Armoricain) délimités par des discontinuités tectoniques, et 
correspondant chacun à des conditions orogéniques ou métamorphiques différentes. 
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Figure 1.1 : Carte simplifiée du Massif Armoricain. A bassin versant de la Vilaine; 
B bassin versant du Blavet. 



2.1.1.1. Domaine Nord-Armoricain 

D'âge essentiellement précambrien et fortement marquée par l'orogénèse cadomienne, 
cette zone est surtout composée de granités fini-cadomiens intrusif s dans des formations 
sédimentaires et volcaniques briovériennes. L'influence de l'orogénèse hercyniennne 
(déformation et mise en place des granités hercyniens) n'apparaît qu'à la périphérie du 
domaine. 

2. 1.1. 2. Domaine Sud-Armoricain 

De structure complexe, cette région est fortement marquée par l'orogénèse 
hercynienne. Elle est essentiellement composée des roches métamorphiques et de 
granitoïdes. Le cisaillement Sud- Armoricain (CSA) grande structure dextre la limite au 
nord. 

2.1.1.3. Domaine Nord-Ouest Armoricain 

Ce domaine, comme le précédent, est marqué par l'orogénèse hercynienne et présente 
principalement des roches métamorphiques et des granitoïdes. 

2.1.1.4. Domaine Médio-Armoricain 

Ce domaine présente, dans sa partie orientale, du Briovérien peu impliqué dans 
l'orogénèse cadomienne et du Paléozoïque discordants (ordoviciens à carbonifères), et 
dans sa partie occidentale, du Briovérien, du Paléozoïque et de nombreux granités 
hercyniens. Les terrains briovériens et paléozoïques sont affectés par des déformations 
d'intensité moyenne contemporaines de la mise en place des granités carbonifères. 

Les bassins versants de la Vilaine et du Blavet sont inclus dans le domaine Médio- 
Armoricain et couvrent approximativement la moitié de ce domaine. Une description très 
générale du substrat géologique de ces deux bassins suffit à les différencier d'un point de 
vue géologique. 

2.1.2. Bassin versant de la Vilaine 

Le bassin versant de la Vilaine est de loin le plus étendu des bassins versants bretons 
(12 000 km 2 ). Les limites septentrionale et méridionale correspondent assez bien aux 
anomalies tectoniques des cisaillements Nord et Sud Armoricain. Les lignes de crête 
occidentale et orientale ne semblent pas soumises à des contraintes de nature lithologique. 
On peut les définir par une ligne géographique passant à l'Est par Laval-Chateaubriand- 
Blain-La Roche Bernard, et l'Ouest par une ligne passant par Qumtin-Rohan-Locrniné- 
Questembert. L'essentiel du bassin versant est caractérisé par du schiste briovérien 
(bassin de Rennes, bassin de Messac, plateau de Rohan) et des sédiments paléozoïques 
(synclinaux du sud de Rennes). Les formations granitiques, localisées au Nord (Quintin, 
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Moncontour) et au sud (Lizio, Questembert), sont niinoritaires et ne représentent qu'une 
faible superficie par rapport à celle du bassin. La zone paléozoïque, constituée de grès 
armoricains, s'étire d'Est en Ouest suivant un axe Nord 1 10° de la forêt de Paimpont aux 
limites orientales du bassin. On retrouve également des formations paléozoïques sur toute 
la partie sud-est du bassin jusqu'au nord de la faille de Nort-sur-Erdre. Les reliefs 
dominants correspondent à des formations granitiques (granités de Quintin, de 
Moncontour, de Lizio) ou paléozoïques (plateau de Guichen). 

2.1.3. Bassin versant du Blavet 

D'une superficie plus modeste (2 000 km 2 ) que celui de la Vilaine, le bassin versant 
du Blavet est néanmoins, de part sa couverture géographique, l'un des principaux bassins 
versants bretons. Accolé à l'Est au bassin versant de la Vilaine, il est limité au Nord par le 
granité de Quintin et à l'Ouest par une ligne géographique passant par Quintin-Guéméné- 
Plouay. Ce bassin présente, comme celui de la Vilaine, une géologie contrastée avec 
cependant un plus grand nombre de massifs granitiques. On y distingue du Nord vers le 
Sud: 

- le granité de Quintin; 

- les séries paléozoïques du synclinal de Chateaulin; 

- le granité de Rostrenen; 

- le granité de Pontivy; 

- enfin à l'Est, les schistes briovériens. 

Les points topographiques dominants sont à nouveau situés sur des formations 
granitiques tels que les granités de Quintin et de Rostrenen. 

2.2. Contexte géomorphologique 

Le relief peu prononcé du Massif Armoricain est relativement bien corrélé avec les 
formations lithologiques sous-jacentes. L'érosion différentielle suivant la dureté du 
substrat a joué un rôle prépondérant dans l'évolution du relief de la péninsule bretonne 
(reliefs appalachiens au Sud du bassin de Rennes). Cependant l'existence de sédiments 
marins récents à diverses altitudes (De Martonne, 1906) met en évidence des anomalies 
géomorphologiques qui témoignent d'événements tectoniques tertiaires et quaternaires 
pouvant être en parie à l'origine du relief actuel. L'histoire du relief breton à partir du 
tertiaire est loin d'être analogue à celle d'une pénéplaine ayant suivi un cycle continu. Les 
grandes étapes de son évolution ont intéressé très tôt les géologues et les géographes (De 
Martonne, 1906; Meynier, 1940; Musset, 1917, 1922, 1928; Guilcher, 1939, 1948, 
1951; Wyns, 1991). On distingue généralement quatre grandes périodes érosives - pré- 
liasique, pré-cénomanienne, paléogène et néogène - conduisant à la formation de surfaces 
pénéplainées et à la mise à nu de structures hercyniennes et précambriennes. 
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La plupart des auteurs s'accordent pour faire débuter la formation du relief actuel au 
Paléogène (De Martonne, 1906; Musset, 1917; Meynier, 1940; Guilcher, 1948; Wyns, 
1991). A l'Eocène, la position élevée de la surface paléogène (+220 à +250 m par rapport 
au référentiel actuel), située au-dessus du niveau marin (+200 à +220 m), autorise le 
nivellement de l'ensemble du Massif Armoricain. A cette époque une pénéplaine met en 
valeur des reliefs comme le Mont d'Arrée, la Montagne Noire et le Menez Hom. Puis 
toujours au Tertiaire, un bombement d'ensemble engendre l'apparition de grabens 
(bassins de Rennes et de Nantes) et de failles formant les reliefs de Hédé, de Pont-Réan 
et du sillon de Bretagne. Suit une succession de phases transgressives intercalées de 
phases plus ou moins érosives, dont les plus importantes se situent au Pliocène et au 
Monastirien. Au Quaternaire, l'évolution péri-glacière marque le relief breton sous la 
forme d'épaisses coulées de solifluxion encore en place de nos jours. L'érosion dans les 
pays tempérés étant relativement faible, on peut considérer que le relief actuel est à peu 
près exactement celui obtenu à la fin de la dernière grande glaciation. 

2.3. Description qualitative des réseaux hydrographiques 

2.3.1. Système hydrographique de la Bretagne 

La lecture de la carte de l'ensemble du réseau hydrographique de la Bretagne révèle 
une organisation générale asymétrique (Fig. 1.2). En effet, l'ensemble des fleuves 
débouchant en Manche voient leur expansion limitée au nord du cisaillement Nord 
Armoricain. De même, en dehors du bassin versant de l'Aulne qui trouve son exutoire 
dans la rade de Brest, la quasi totalité des domaines Médio et Sud Armoricain est drainée 
par des fleuves débouchant sur la façade atlantique. Il en résulte une pente globale du 
domaine Nord-Armoricain du Sud vers le Nord, et des domaines Médio et Sud 
Armoricain du Nord vers le Sud. Un autre trait marquant du réseau breton est 
l'alignement très net des cours d'eau suivant le cisaillement sud-armoricain. Les 
orientations des rivières et des fleuves, au Nord et au Sud du cisaillement, montrent 
cependant des sens d'écoulement initiaux Nord-Sud et non pas suivant une orientation 
Est-Ouest, comme le laissent suggérer certaines configurations actuelles. Nous observons 
donc une réorientation des sens de drainage postérieure à l'installation de l'ensemble du 
système hydrographique breton. Ceci laisse supposer un mouvement tectonique récent du 
domaine sud armoricain, qui se matérialise au niveau des zones de faiblesse constituées 
par le cisaillement Sud armoricaine. 

Si les structures géologiques ont une influence certaine sur le cours de certaines 
rivières, il est néanmoins très surprenant de constater la constance des cours d'eau 
principaux à suivre des directions d'écoulement pour le moins paradoxales. L'axe général 
d'écoulement des rivières principales, l'Hle-Vilaine, la Rance, le Couesnon inférieur, le 
Blavet, l'Oust, l'Oudon, la Mayenne, suit la direction des méridiens. Or dans bien des 
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cas, les structures topographiques s'étirent suivant une direction Est-Ouest: crête de 
Bécherel-Hédé-St Aubin du Cormier, bassin de Rennes, synclinaux paléozoïques du Sud 
de Rennes. On note donc une discordance entre certaines caractéristiques du réseau 
hydrographique et le système orographique existant. Ceci peut s'expliquer simplement 
par des phénomènes de surimposition (Fig. 1.3; Meynier, 1940) des axes principaux 
d'écoulement, antérieurs aux reliefs actuels. 

2.3.2. La Vilaine et ses affluents 

La Vilaine, cours d'eau majeur du Massif Armoricain, d'une longueur propre de 220 
km mais représentant un réseau hydrographique de plus de 4000 km, n'échappe pas aux 
remarques précédentes. 

Au Sud de Rennes, l'alternance de roches dures (grès armoricains) et de roches plus 
tendres (schistes briovériens et paléozoïques) vis à vis de l'érosion, (synclinaux 
paléozoïques) est à l'origine de l'émergence d'un relief appalachien au centre du bassin 
versant de la Vilaine (De Martonne, 1906). Le réseau hydrographique associé suit bien 
évidemment les mêmes tendances et présente un alignement des structures de drainage (le 
Meu, la Seiche, le Semnon, le Canut). L'orientation Est-Ouest des plis affectant les 
terrains paléozoïques, perpendiculaire à l'orientation globale du bassin versant, provoque 
également l'apparition d'une véritable entité hydrographique au sein du réseau 
hydrographique de la Vilaine, sous la forme d'une structure étoilée drainant le bassin de 
Rennes. Pourtant ces mêmes plateaux ne peuvent résister au pouvoir érosif de la Vilaine, 
qui les traverse de part en part par un processus de surimposition (Fig. 1.4; De Martonne, 
1906), et révèlent ainsi une configuration antérieure du réseau hydrographique bien 
différente de l'actuelle. Le tracé du cours des affluents de la Vilaine et la position d'un 
certain nombre de cluses dans les synclinaux paléozoïques du Sud de Rennes permettent 
de reconstituer l'ancien réseau hydrographique (Fig. 1.5; Meynier, 1940). Ce dernier se 
développait dans le sens Nord-Sud, avant d'être peu à peu déformé par un jeu de captures 
favorisant les drains des roches les plus tendres. L'étude de la Vilaine supérieure et de la 
basse Vilaine révèle une évolution similaire du réseau hydrographique (Meynier, 1940). 

La Vilaine est composée schématiquement de quatre réseaux successifs de formes 
différentes: la Vilaine supérieure (réseau en épi dans des formations briovériennes), les 
synclinaux du Sud de Rennes (réseau en espalier dans des formations paléozoïques), la 
Vilaine du bassin de Messac (réseau en épi dans des formations briovériennes), et la 
Vilaine méridionale (réseau en espalier dans des formations granitiques et paléozoïques). 
Cette dichotomie du réseau reflète le contrôle progressif des cours d'eau par la géologie 
sous-jacente qui déplace les confluents vers l'amont. Mais la persistance de certains traits 
originels (cours de la Vilaine, de l'Ille et de Isac) témoigne également d'une relative 
indépendance du réseau vis à vis de son environnement actuel. 
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Fig l.3:Processus de suri mposition. 




Fig l.4:Diagramme des plateaux siluriens 
traversés par la Vilaine (De Martonne, 1 906) 



Vilaine 

ancien cours d'eau 
-~^^y cours d'eau actuel 

10 km 

■ i - _ limite sud du bassin de Rennes 



1 C 



cluse, seuil 



Fig 1.5 : Modification du cours d'eau de la Vilaine et de 
ses principaux affluents due à l'émergence d'un relief 
appalachien au Sud du Bassin de Rennes 
(d'après Meynier, 1940). 



2.3.3. Le Blavet 

Le réseau hydrographique du Blavet, moins étendu, présente moins de structures 
anecdo tiques. On peut simplement évoquer l'influence des granités de Pontivy et de 
Rostrenen qui tendent à dévier le cours du Blavet vers l'Est. 

2.4. Contexte climatique 

Il existe peu de monographies sur les caractéristiques hydro-climatiques de la 
Bretagne. On peut néanmoins citer une référence de 1969 de Larivière et Verdou 
résumant les principaux traits du climat breton des années 50 à 70. 

La péninsule bretonne de part sa situation géographique bénéficie de l'influence 
océanique des courants marins (Gulf-Stream) et aériens (vents d'ouest associés à la 
dépression d'Islande). Son climat typiquement océanique se caractérise par des hivers 
cléments, des étés frais, et de faibles amplitudes de température. Les précipitations sont 
inégalement réparties sur le territoire (Fig. 1.6). La quantité de pluie annuelle est 
globalement croissante d'Est en Ouest: 650 mm à Rennes, 800 mm à Pontivy et 1000 mm 
dans le Finistère. On note cependant une forte influence du relief local: 1500 mm sur les 
reliefs finistériens et 800 mm dans le bassin de Pontivy situé à l'abri de la Montagne 
Noire. Le littoral notamment Sud subit également un régime particulier avec une 
pluviométrie annuelle relativement faible de 700 mm: les nuages se formant quelques 
kilomètres à l'intérieur des terres. La carte du déficit en eau (Pluie-ETP ) pour les 
périodes d'été (Fig. 1.7) montre un fort déficit de réserve en eau dans la partie orientale de 
la Bretagne. 

Les bassins versants du Blavet et de la Vilaine sont soumis à des régimes climatiques 
différents tant sur le plan pluviométrique que sur le plan du déficit en eau. Le bassin 
versant du Blavet présente en effet un déficit en eau («-150 mm) qui est globalement 
inférieure à celui de la Vilaine (=-250 mm). 

2.5. Utilité de l'analyse quantitative 

L'exemple presque caricatural de configuration hydrographique de la Vilaine, a le 
mérite de mettre en évidence certaines particularités de la dynamique des réseaux 
hydrographiques, et par là-même des processus d'érosion. 

Un réseau hydrographique est le produit complexe d'une succession de conditions 
limites et locales qui peuvent changer radicalement son évolution. Cependant, certains 
cours d'eau acquièrent, de part leur histoire, une puissance érosive qui leur donne une 
indépendance vis à vis des conditions locales. Le rapport de l'amplitude des perturbations 
sur les conditions locales ou aux limites (changement d'induration des roches, variations 



16 



Septembre - Octobre - Novembre 




Décembre - Janvier - Février 





Mars - Avril - Mai 



Fig 1.6 : Hauteur moyenne des précipitations (période d'étude 1951 - 1975) 
sources: MétéoFrance 




Fig 1.7 : Bilan (Pluie - Evapotranspiration Potentielle) sur la période du 1/3 au 31/8 
(période d'étude 1951 - 1975); Les zones les plus sombres sont celles pour 
lesquelles le risque de sécheresse (au printemps été) est important. 
(Raphalen - I.T.C.F, 1982) 



Evapotranspiration Potentielle (ETP) : lorsque le seul facteur est le climat 
l'ETP est définie comme étant la quantité d'eau maximale susceptible de 
s'évaporer d'un gazon ras recouvrant totalement le sol, bien alimenté en eau, en 
phase active de croissance et situé au sein d'une parcelle étendue. En conditions 
naturelles, aux facteurs climatiques peuvent s'ajouter d'autres facteurs limitant 
l'évapotranspiration (réserve en eau du sol, état de développement des 
végétaux). On parle alors d'avapotranspiration réelle (ATER). La fraction d'eau 
retournée à l'atmosphère varie entre 50% pour les régions tempérées et 
pratiquement 100% en régions arides (source : Bruneau, 1992). 



climatiques, actions anthropiques etc...) et du pouvoir érosif du cours d'eau, détermine le 
devenir de l'orientation de ce dernier et par là-même de la distribution spatiale des aires 
drainées. Les phénomènes de surimposition par exemple marquent une non linéarité des 
processus et laissent suggérer l'existence de seuils critiques de comportement. Au-delà 
d'une valeur seuil de taux d'érosion, le réseau est indifférent à toute perturbation et 
s'organise suivant des règles physiques encore non expliquées. En-deçà du seuil, 
l'instabilité peut engendrer des modifications plus ou moins radicales sur les conditions 
d'écoulement et marquer définitivement l'évolution du cours d'eau (capture par 
déversement, Fig. 1.8). La morphologie d'un bassin versant, étroitement liée à 
l'extension du réseau hydrographique qui le draine, suit le même type d'évolution. Un 
bassin versant, sur une échelle de temps long, tend vers un état d'équilibre (pénéplaine). 
Ce dernier sans cesse perturbé par des variations de nature tectonique, climatique ou 
géomorphologique (extension des bassins versants voisins) voit son extension remise en 
question, suivant l'amplitude des perturbations et sa surface de couverture géographique. 
Comprendre l'origine de la stabilité ou de l'instabilité des structures géomorphologiques 
nécessite au préalable la compréhension des phénomènes physiques contrôlant l'évolution 
de tels systèmes. 

La géométrie de la topographie et des réseaux hydrographiques traduit, par des 
structures visibles, la nature des phénomènes mis en jeu. Trouver les lois d'organisation 
géométrique qui régissent ces systèmes constitue une première étape dans la 
compréhension globale de la physique qui leur est rattachée. L'approche quantitative 
permet d'établir de manière plus formelle les règles d'organisation de ces systèmes 
naturels, et par là-même d'aborder le problème de leur modélisation. 

3. Codification et nomenclature 

3.1. Définitions 

L'analyse de la structure des réseaux hydrographiques emploie un certain nombre de 
définitions, afin de caractériser les éléments de rivière qui vont constituer les populations 
statistiques. La terminologie utilisée dans ce travail s'inspire fortement de celle établie par 
les auteurs anglo-saxons, précurseurs en ce domaine (Horton, 1945; Strahler, 1964; 
Shreve, 1966). Afin de permettre au lecteur de se familiariser avec la terminologie utilisée 
dans les revues internationales, les termes spécifiques seront suivis entre parenthèses de 
leurs équivalents en anglais. 

Une arborescence hydrographique est composée de portions de cours d'eau ou biefs. 
Ce terme très général ne permet pas à lui seul de définir les différents éléments du 
système. Aussi doit-on spécifier plus précisément la nature du bief en définissant ses 
limites géographiques. La démarche la plus naturelle est de considérer les biefs délimités 
par des confluences (fork) ou des sources (source). Shreve, dans son introduction au 
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capture par déversement 




capture par recul de tête 



Fig 1.8 : Processus de capture. 



Random Topology Model (1966), propose une classification simple et indépendante de 
toute règle artificielle. Les biefs y sont classés en deux groupes: les biefs intérieurs 
{interior links) reliant deux confluences successives dans l'arborescence de la rivière, et 
les biefs extérieurs (exterior links, nombre de sources drainant en un point donné du 
réseau: magnitude ) reliant les sources aux premières confluences d'aval (Fig. 1.9). Cette 
distinction trouve son origine dans les résultats observés dans les années 50 (Strahler, 
1954; Schumm, 1956; Melton, 1957; Morisawa, 1962), où les biefs extérieurs et 
intérieurs semblent suivre des distributions statistiques différentes. 

3.2. Codification de Strahler 

Un observateur étudiant la structure du réseau hydrographique ressent intuitivement le 
besoin de classer, suivant une certaine hiérarchie, les différents biefs composant 
l'arborescence de la rivière. Ceci afin d'attribuer, dans la description de chaque segment 
de cours d'eau, un poids relatif au débit qui le traverse. Un certain nombre de 
codifications sont proposées dans la littérature et présentent un intérêt plus ou moins 
grand suivant l'objectif fixé (codification de Horton, de Strahler, méthode de "priorité à 
gauche"). 

Pour analyser le réseau hydrographique, nous avons utilisé la codification introduite 
par Strahler en 1952. Le choix de la codification de Strahler pour ce travail, vient 
principalement du souci de pouvoir se référer aux travaux antérieurs, qui dans leur grande 
majorité utilisent cette codification. Strahler définit l'ordre d'un bief de la façon suivante 
(Fig. 1.9) (bief défini suivant la codification de Strahler: stream): 

• tout bief sans affluent est d'ordre 1; 

• un bief formé par la confluence de deux biefs d'ordre (o est d'ordre ctw-l; 

• un bief formé par la confluence de deux biefs d'ordres différents prend l'ordre du 
bief le plus élevé: co +co' -> max(ô),û/). 

On définit un réseau hydrographique par la valeur de l'ordre le plus élevé de son 
arborescence. Avec le réseau hydrographique vectorisé au 1/100 000, la Vilaine aurait 
ainsi un ordre de 7 et le Blavet, un ordre de 6. 

Si la méthode de Strahler permet de classer hiérarchiquement les biefs en fonction de 
leur degré de connexion, elle n'en reste pas moins très dépendante de la précision avec 
laquelle les biefs sources sont définis. De plus, la règle (û) +Cû' -> max(û), û)J) n'est pas 
associative, et ne traduit pas l'augmentation de débit induite par l'apport en eau du bief 
d'ordre le plus faible. La hiérarchie suivant la codification de Strahler ne rend donc pas 
compte des processus d'intégration sur des variables géométriques, comme l'aire drainée, 
ou physiques, comme le débit. Elle impose seulement une organisation dans l'ensemble 



21 



ext 



int 



ext: biefs extérieurs 
int: biefs intérieurs 



int 



,ext 



Al, A2: aire drainée par un bief 
A3, A4: aire drainée propre d'un bief 




Codification de Strahler 



Fig 1.9 : Nomenclature de biefs 



des connexions d'un réseau hydrographique. Aussi est-il hasardeux d'établir un parallèle 
entre une hiérarchie des connexions et une hiérarchie des débits. 

4. Traitements des données 

4.1. Réseau hydrographique vectorisé 

4.1.1. Source cartographique 

Les réseaux hydrographiques de la Vilaine et du Blavet ont été vectorisés à partir des 
cartes thématiques de l'Agence de bassin Loire-Bretagne, elles-mêmes obtenues sur la 
base des cartes IGN au 1/100 000 datées de 1986. Nous avons retenu comme 
information, le réseau représenté par un trait plein continu et éliminé tant que possible tout 
drain d'origine anthropique. 

La vectorisation du réseau cartographié a l'avantage de fournir une donnée issue 
directement d'observations terrain et de photographies aériennes, et indépendante de tout 
modèle de construction de réseau de drainage à partir d'une représentation numérisée de 
la topographie (voir § suivant). Elle présente cependant l'inconvénient d'être imprécise 
dans la détermination des biefs sources. En effet, la position géographique des têtes de 
biefs varie fortement suivant l'état hydrique des bassins versants d'amont. Un bief 
d'ordre 1 peut voir sa longueur multipliée par 2 ou 3 et devenir un bief d'ordre 2, au 
cours d'un orage. De plus, dans bien des cas les sources ne sont pas délimitées par un 
point géographique précis, mais se présentent sous la forme de zones saturées en eau 
réparties de façon plus ou moins discontinue le long de la ligne de thalweg (Montgomery 
et Dietrich, 1988). Le report de la position des sources est donc très dépendant du 
moment auquel le géographe effectue le relevé et de son appréciation. La position 
géographique des biefs sources dans le fichier vectorisé est par conséquent sujet à 
caution. 

4.1.2. Vectorisation 

La vectorisation est réalisée par le logiciel Descartes conçu à L'atelier de Spatialisation 
Numérique de l'Ecole Nationale Supérieure d'Agronomie de Rennes. Ce logiciel fournit à 
partir d'un fichier raster, un fichier vecteur où chaque bief délimité par deux confluences 
ou par une source est individualisé par un fichier de points qui lui est propre. Les biefs 
peuvent ainsi être manipulés indépendamment des uns des autres. Si l'on tient compte des 
déformations dues au support papier d'origine et à la retranscription sur un document 
propre, l'ensemble de la vectorisation permet d'obtenir un réseau vectorisé dont la 
précision est de l'ordre de 100 à 200 mètres (Fig. 1.10). 
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Fig. 1-10 : Réseaux hydrographiques vectorisés au 1/100 000. 

Abscisse : longitudes en Lambert 2 étendu. 

Ordonnées : latitudes en Lambert 2. a) Vilaine, b) Blavet. 



4.1.3. Construction d'un réseau simple 

Le réseau naturel présente, notamment dans les zones de faible pente, des îles et des 
bifurcations dans le sens amont-aval. Cette caractéristique anodine des rivières, n'en est 
pas moins embarrassante pour le numéricien qui préfère travailler avec des arborescences 
parfaites, à savoir sans bifurcation amont-aval (Fig. 1.11). La plupart des codifications et 
des modèles de réseaux de drainage assimile implicitement les réseaux hydrographiques à 
des arborescences parfaites (ex: codification de Strahler). Si cette représentation 
correspond effectivement à la réalité pour une certaine échelle, elle reste beaucoup trop 
simpliste pour reproduire le tracé du lit mineur d'un cours d'eau. Cependant, pour 
pouvoir appliquer la codification de Strahler et simplifier le traitement numérique du 
réseau vectorisé, nous avons préféré extraire une arborescence parfaite du réseau initial. 

Un module de calcul conçu à cette occasion, permet d'effectuer ce travail de façon 
automatique, quel que soit le réseau vectorisé. La procédure consiste dans les grandes 
lignes à déterminer les biefs formant des mailles fermées, et à éliminer l'un des biefs de la 
maille. Le choix du bief à éliminer devrait tenir compte du débit qui lui est associé. Le 
fichier vectorisé ne comporte que des informations sur la géométrie du bief et ne permet 
donc pas de choisir un critère pertinent de sélection. Aussi nous avons préféré garder une 
procédure aléatoire (ordre d'apparition dans la filiation des biefs) pour le choix du bief à 
retirer. 

La codification de l'arborescence suivant les règles de Strahler fait aussi l'objet d'un 
traitement numérique. Un autre module de calcul définit, à partir des biefs de base 
(portion de cours d'eau entre deux confluences), les biefs et l'ordre auquel ils 
correspondent suivant cette codification. 

On peut dès lors, à l'issue de ces deux traitements, effectuer la statistique de la 
géométrie des biefs appartenant à une classe donnée. 

4.2. Traitement numérique de la topographie digitalisée 

Avant l'avènement des moyens informatiques actuels, le calcul de l'aire drainée 
{drainage area), pour des raisons de lourdeur de calcul, se limitait à des bassins versants 
inférieurs à 100 km 2 (Strahler, 1954; Schumm, 1956; Melton, 1957; Morisawa, 1962). 
L'apparition de topographies digitalisées {DEM: Digital Elévation Model) a permis de 
concevoir des logiciels de calcul automatique d'aire drainée et d'accéder ainsi au 
traitement de plus grandes surfaces. 
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réseau vectorisé d'origine réseau vectorisé "nettoyé" 

Fig 1.11 : Transformation du réseau vectorisé d'origine en une arborescence simple 
(suppression des confluences aval est) 




Fig 1.12 : Limite géographique d'un bassin versant (d'après Roche 1971) 



4.2.1. Modèle Numérique de Terrain 
On distingue deux types de MNT: 

- les MNT à mailles rectangulaires. Ds ont le principal avantage d'être d'utilisation 
assez simple (une seule valeur par point) et facilement superposables aux images 
satellitaires. Ils présentent l'inconvénient de n'être pas adaptés à des applications 
hydrologiques (mauvaise détermination de la pente locale réelle). On peut noter par 
exemple, l'existence de minima locaux artificiels (pits ou sinks) nécessitant un pré- 
traitement du fichier assez lourd. 

- les MNT à mailles triangulaires (TIN: Triangular Irregular Network). 
L'échantillonnage des altitudes s'effectue suivant la convexité des courbes de niveau et 
n'est par conséquent pas régulier. Cette procédure retranscrit assez fidèlement la 
distribution des pentes locales, mais l'utilisation de tels fichiers d'altitudes peut s'avérer 
fastidieuse (Palmer, 1984; Frank et al., 1986). 

La suite du texte fera implicitement référence à des MNT à mailles rectangulaires, ces 
derniers étant les plus couramment utilisés dans les applications hydrologiques. La 
détermination de l'altitude, en chaque noeud de la grille d'échantillonnage, est obtenue à 
l'aide de modèles d'interpolation de cotes de référence (courbes de niveau). Ces modèles 
peuvent difficilement s'adapter à toutes les configurations topographiques et présentent 
des erreurs d'estimation plus ou moins fortes suivant les cas de figure. L'erreur maximale 
est égale au demi intervalle entre chaque courbe de niveau du document cartographique 
d'origine (si la différence de dénivelé entre deux courbes de niveau est de 5 mètres, 
l'erreur maximale sur les altitudes sera de 2,5 mètres). Des méthodes 
photogrammétriques appliquées à des vues aériennes ou des images SPOT offrent une 
plus grande densité de points de référence et permettant d'affiner restimation. Mais elles 
exigent une certaine prudence dans l'interprétation des altitudes, car le signal brut utilisé 
mesure l'altitude du toit du couvert végétal et des bâtiments. 

4.2.2. Notion d'aire drainée 

On appelle aire drainée en un point l'aire délimitée par le périmètre à l'intérieur duquel 
toute l'eau précipitée se dirige vers ce point. Lorsque ce point appartient au cours d'eau, 
l'aire drainée représente alors le bassin versant (watershed, en américain et catchment en 
anglais) du point d'étude. Cette définition quoique très simple est très dépendante du 
système hydrologique que l'on considère. En effet, pour un même point, l'aire drainée 
peut être très différente si l'on tient compte ou non des écoulements souterrains (Fig. 
1.12). Les limites géographiques de l'aire drainée sont le plus souvent assimilées au 
secteur défini par les deux lignes de plus grande pente aboutissant au point exutoire 
(outlef) et par la ligne de crête qui les rejoint. Cette définition exacte dans le cas d'eau 
ruissellée, ne tient pas compte de l'étendue géographique réelle des échanges entre les 
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différents compartiments composant l'hydrosphère (atmosphère, surface, sol et sous- 
sol). Mais pour des raisons de simplification, de disponibilité de l'information, de 
logique et de précision suffisante (variable suivant l'échelle du bassin), les hydrologues 
utilisent généralement cette définition. 

4.2.3. Calcul de l'aire drainée 

Le calcul de l'aire drainée en un point donné du système passe par la détermination de 
la direction d'écoulement en chaque maille du MNT. Une fois cette opération accomplie, 
il est possible de définir le réseau de drainage constitué par l'ensemble des chemins 
d'écoulement. Les logiciels de construction de réseau de drainage se différencient 
principalement suivant le mode de répartition des écoulements au niveau de chaque 
maille. La méthode la plus couramment utilisée est la méthode D8 introduite par 
O'Callaghan et Mark (1984). Dans cette méthode chaque pixel se déverse entièrement 
dans l'un de ses 8 voisins suivant le critère de plus grande pente (Fig. 1.13). Cette 
procédure, simple à mettre en oeuvre, reste très rigide dans le choix des directions de 
drainage possibles (8 directions possibles) et introduit un biais important dans 
l'estimation de l'aire drainée pour des directions intermédiaires d'écoulement. Afin de 
minimiser cette erreur systématique, Fairfield et Lemayrie (1991) proposent une nouvelle 
méthode (Rh8) en multipliant la pente locale par une composante stochastique. La 
méthode de Léa (1992) basée sur la pente du plan passant au mieux par les 4 noeuds de 
chaque maille, offre également la possibilité de mieux définir les directions de drainage, et 
présente l'avantage, par rapport à la méthode Rh8, d'être reproductible pour un état 
donné du système. 

Les modèles précédents attribuent une seule direction de drainage à chaque pixel et 
conduisent à la création de réseaux arborescents sans bifurcation amont-aval. Si cette 
représentation correspond aux systèmes naturels à une échelle globale, elle n'en reste pas 
moins peu crédible pour certaines configurations topographiques. Les systèmes naturels 
présentent en effet des structures qui témoignent d'une distribution de volumes d'eau 
multidirectionnelle (réseaux en tresse, delta, cônes de déjection). Pour compenser cette 
lacune, Quinn et al. (1991) d'une part, et Freemann (1991) d'autre part, proposent une 
distribution du débit pondérée par une fonction de la pente locale, dans les pixels voisins 
d'altitude inférieure. La répartition systématique des débits engendre cependant une 
exagération de l'étendue des aires drainées. Enfin Costa-Cabral et Burges (1994) 
présentent dans leur modèle DEMON, une synthèse des différents points forts des 
modèles précédents avec une représentation en deux dimensions des écoulements au 
niveau de chaque pixel. 

Nous avons utilisé pour ce travail la méthode D8. Les inconvénients liés à l'existence 
d'une direction unique sont ici minimisés par la taille relativement importante du pixel 
(250 x 250 m). On peut en effet considérer que le critère de plus grande pente est une 
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Fig 1.13: Modèle de drainage et calcul de l'aire drainée 
(critère de plus grande pente). 
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Fig 1.14: Correction des minima locaux 
(critère de plus grande pente) 



estimation raisonnable de la direction moyenne d'écoulement, pour des pixels 
représentant 25 à 50 % de la surface moyenne d'un versant (cas du MNT de la Bretagne 
au pas de 250 m). 

4.2.4. Les minima 

Suivant la qualité et le pas du MNT, l'application des règles précédentes de drainage 
conduit à l'apparition d'un pourcentage d'exutoires non raccordés à la mer (points dont 
tous les voisins ont une altitude supérieure ; Fig. 1.14). Si certains de ces minima ont une 
réalité physique (dolines des karsts ou lacs naturels), ils sont dans leur grande majorité 
des artefacts liés à la mauvaise représentativité de l'échantillonnage des altitudes. La 
plupart des modèles de construction de réseau de drainage part du principe que tous les 
minima sont artificiels, et utilise des techniques spécifiques pour les supprimer. La 
technique la plus courante est de rechercher sur le périmètre du bassin versant artificiel, 
créé par le minimum à supprimer, le col le plus bas et d'orienter les directions de drainage 
vers ce col (Mark, 1983; Mark et al., 1984; Jenson et Domingue, 1988; Morris et 
Heerdegen, 1988). 

L'ensemble des opérations numériques réalisant le calcul de l'aire drainée est effectué 
par le logiciel MNTSurf conçu à L'atelier de Spatialisation Numérique de l'Ecole 
Nationale Supérieure d'Agronomie de Rennes. 

4.2.5. Utilisation du réseau vectorisé dans le calcul de l'aire drainée 

Nous avons modifié le programme d'origine, afin de pouvoir superposer au MNT le 
réseau vectorisé dont nous disposons, et améliorer ainsi la précision du réseau de 
drainage résultant 

La méthode que nous avons suivie consiste dans un premier temps à définir les pixels 
du MNT comprenant au moins un point du réseau vectorisé. On définit ainsi un ensemble 
de "pixels réseau" pour lesquels la direction de drainage est imposée par le fichier 
vectorisé. Le MNT présente, suivant sa précision (pas de la maille d'échantillonnage), 
une erreur plus ou moins forte dans la détermination de la position géographique réelle 
des talwegs. Cette imprécision engendre le positionnement de "pixels réseau" sur les 
versants et non pas dans le fond des vallées. De manière à adapter l'information 
vectorisée à la précision du MNT, nous contrôlons pour chaque "pixel réseau" son 
altitude par rapport à celles des pixels voisins non "pixels réseau" et mitoyens aux "pixels 
réseau" amont et aval (Fig. 1.15). Si l'un des voisins possède une altitude inférieure, le 
"pixel réseau" est déplacé sur ce dernier et les directions de drainage modifiées. 
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Superposition du réseau vectorisé 



Correction des pixels réseaux 
en fonction de la topographie numérisée 
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Fig. 1.15 : Superposition du réseau hydrographique vectorisé sur la topographie digitalisée. 
La définition de la topographie digitalisée et du réseau vectorisé étant différentes (précision 
du MNT= 250 m, précision du réseau vectorisé 40 m), certaines configurations de drainage 
imposées par le réseau sont incompatibles avec la grille d'altitude. Une phase de recalage 
des pixels "réseau" est donc nécessaire. 



La prise en compte de la nature des pixels ("réseau" ou "non réseau") permet 
d'améliorer la procédure d'élimination des minima locaux. Après la détermination des 
directions de drainage brutes (sans correction des minima locaux), tous les minima locaux 
sont répertoriés. Avant d'effectuer la correction de ces points aberrants par la méthode 
utilisée par le logiciel, une routine détermine pour chaque minimum local si dans son 
voisinage immédiat, il n'existe pas un ou des "pixels réseau". Si tel est le cas, on 
affectera au minimum une direction de drainage de manière à le raccorder au "pixel 
réseau" voisin d'altitude la plus basse. Cette procédure permet de réduire le temps de 
calcul nécessaire au traitement de l'ensemble du MNT, par une réduction de 30 à 40 % du 
nombre initial de minima. 

5. Méthodes analytiques 

L'analyse quantitative des réseaux hydrographiques s'appuie pour l'essentiel sur 
l'étude des distributions des populations de paramètres topographiques Oongueurs, aires 
drainées) et les lois empiriques reliant ces paramètres entre eux. Avec le concept fractal 
(Mandelbrot, 1975), les géomorphologues disposent d'un outil supplémentaire pour 
appréhender la géométrie des formes naturelles. Nous développerons dans ce chapitre les 
méthodes de calcul liées à ces deux types d'analyse. 

5.1. Calcul des distributions de populations 

5.1.1. Définitions 

Les lois empiriques utilisées en géomorphologie quantitative, telles que les lois de 
Horton, ne se réfèrent qu'à la moyenne des paramètres. En fait, la moyenne peut occulter 
des distributions beaucoup plus complexes comme le montrent, par exemple, les figures 
2 et 3 de l'article constituant la dernière partie de ce chapitre. Nous avons donc étudié 
spécifiquement les distributions de longueurs de drains et d'aires drainées en fonction des 
différents types de biefs: intérieurs, extérieurs et d'ordre co. Nous présentons (i) la 
densité de distribution des longueurs de biefs, n(l,co), où n(l,(o)àl représente le nombre 
de biefs d'ordre œ et de longueur comprise entre / et l+dl, et (ii) la distribution cumulée 
C(l,Cû), égale au nombre de biefs de longueur supérieure à /. C(l,(û) est l'intégrale de / à 
l'infini de la fonction n{l,Cû). 



L'étude comparée de n(l) et C(l) permet en effet d'établir avec plus de précision les 
domaines de validité des distribution théoriques. C(l) donne une bonne estimation de la 
tendance globale de la distribution. La valeur de C(l) est cependant sensible aux effets de 
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taille finie, ce qui peut biaiser l'interprétation. Les variations de n(l) permettent d'une part 
d'estimer l'ordre de grandeur de la précision de la distribution, et d'autre part de mettre en 
évidence les effets de taille finie. 

5.1.2. Calcul et interprétation 

La méthode la plus intuitive pour calculer n(t) consiste à diviser le champ des 
longueurs par un intervalle régulier dl et compter le nombre de biefs pour chaque dl. Le 
choix de la valeur de dl est capital car il détermine le degré de lissage des données 
initiales. La valeur idéale de dl dépend de la précision des données et de la distribution 
théorique escomptée. L'originalité de la méthode utilisée pour ce travail réside dans la 
détermination de la valeur de dl, pour laquelle n(J) présente le moins de fluctuations, et 
ceci pour différentes valeurs du champ des longueurs (Davy, 1993). Certaines 
populations, telles que la longueur totale de réseau à l'amont d'un point, présentent un 
champ très large de longueurs (0.1 - 4000 km). Pour tenir compte des variations 
importantes des domaines des populations, n(J) et C(l) sont calculés pour une variation 
linéaire des longueurs (différence constante entre deux longueurs consécutives) et pour 
une variation logarithmique (rapport constant entre deux longueurs consécutives). 

Une détermination indépendante de C(/) et de n(l) nous permet de définir la loi la plus 
représentative des données (Davy, 1993). Dans le cas de distributions exponentielles 
(exp(-lflo)) ou de lois de puissance (l' a ), le rapport n(l)/C(l) est un moyen puissant pour 
déterminer la nature de la tendance théorique. n(l)/C(l) est en effet égal à 1 pour une loi 
exponentielle et à l' 1 pour une loi de puissance. 

5.2. Formes fractales 

Les formes naturelles telles que la surface des continents ou des océans, les 
arborescences de toute nature, ou encore les nuages, frappent l'observateur par 
rirrégularité d'aspect et ceci à des d'échelles différentes. On peut dès lors s'interroger sur 
la nature des phénomènes physiques mis en jeu et si ces derniers sont, comme les formes 
qu'ils engendrent, similaires quelle que soit l'échelle d'étude. Afin de répondre à ces 
deux questions, on peut rechercher les caractéristiques géométriques intrinsèques de ces 
objets, et établir les modèles physiques possibles engendrant de telles propriétés. La 
géométrie classique étant impuissante à décrire des formes irrégulières quelle que soit 
l'échelle, la création de nouveaux outils mathématiques a dû être nécessaire. 

5.2.1. Définitions et propriétés des fractals 

Une géométrie fractale permet de décrire des objets dont rirrégularité et la complexité 
des formes ne peuvent pas être décrites dans espace géométrique euclidien de dimension 
entière (droite=l, surface=2, volume=3). Mandelbrot (1975) a introduit les dimensions 
fractales pour décrire de tels objets. Prenons l'exemple d'une courbe de géométrie très 
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irrégulière. Plus son degré de sinuosité et d'irrégularité sera important, plus la courbe 
aura tendance à couvrir la surface sur laquelle elle repose. L'espace recouvert peut donc 
avoir une dimension différente de la dimension de la courbe (d=l). Une courbe dense va 
occuper un espace de dimension 2, c'est à dire une portion de plan. Une courbe 
présentant des irrégularités à toute échelle peut occuper un espace de dimension 
intermédiaire entre le plan (d=2) et la ligne (d=l). De plus, si la fonction d'échelle qui 
définit l'amplitude des irrégularités est une loi de puissance, on peut définir une 
dimension fractionnaire pour l'espace recouvert; l'objet est alors appelé fractal. 

On peut exprimer la définition ci-dessus de façon plus formelle à partir de la notion de 
mesure. Mesurer un objet consiste à compter le nombre d'unité de mesure jti (unité de 
longueur, de surface ou de volume) nécessaire au recouvrement total de l'objet. L'unité 
de mesure \l peut être décrite sous la forme générale de sphères de dimension d, ceci 
quelle que soit d. La mesure Me de l'objet E s'exprime alors sous la forme: 



avec N le nombre d'unité de mesure, e l'étalon de mesure (ou le rayon de la sphère 
unité), et d la dimension de l'unité de mesure (d=Q pour un point, d=l pour une droite, 
d=3 pour un volume). La mesure exacte de l'objet est obtenue pour £-»0. Hausdorff et 
Besicovitch ont généralisé la définition d'une mesure en introduisant la notion de sphère 
unité de dimension fractionnaire. Es ont ainsi démontré que pour tout objet, il existe une 
dimension D telle que: 



Dans le cas d'objets simples D est égale à la dimension topologique de l'objet (dimension 
de l'objet en géométrie euclidienne) i.e. 0 pour un point, 1 pour une ligne, 2 pour un 
plan... D peut également prendre des valeurs non-entières pour des objets de géométrie 
irrégulière quelle que soit l'échelle. De tels objet sont appelés fractals. 

5.2.2. Méthodes de calcul 

Pour des raisons pratiques, les méthodes de calcul de la dimension fractale d'objets 
complexes n'utilise pas la notion de sphères unités de dimension fractionnaire. Si l'on 
recouvre un objet fractal avec des sphères de dimension d entière, le nombre de sphères 
nécessaires N(e) dépendra de la valeur e du rayon des sphères. On montre alors que: 



M E = N/i = Né** (1) 



lim £ ^o M E ->°° 



d<D 



lim e ^oM£ -> 0 



d>D 



lim g ^o N(e) ~ e 



-D 



(2) 



D = lim 



Log(N(£)) 
Log(e) 



(3) 
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La loi de puissance sur la mesure définit une certaine auto-similarité de l'objet. En effet, 
le rapport de longueur pour deux longueurs étalons différentes ne dépendra que du 
rapport des longueurs étalon et non de leurs valeurs absolues; autrement dit, quel que soit 
l'agrandissement utilisé pour regarder la courbe, on retrouvera des formes dont la 
géométrie est statistiquement similaire. 

La dimension fractale peut donc être mesurée en suivant l'évolution de N(£) en 
fonction de l'échelle £. La mesure concrète de la dimension fractale consistera à recouvrir 
l'objet de sphères de rayon e, à compter le nombre de sphères, et à réitérer l'opération 
pour un e inférieur. Ces méthodes de recouvrement s'effectuent généralement avec des 
sphères de dimension égale à celle de l'espace englobant. On utilise essentiellement deux 
méthodes de recouvrement: 

le comptage des boîtes (box counting) où les sphères ont une dimension égale à 
celle de l'espace englobant, 

la méthode du compas (compass method) où la dimension des sphères est égale 
à 1 (segment de droite e). 

52.2.1. Méthode des boîtes (exemple de sphères unité de dimension 2) 

L'espace d'étude, représenté sous la forme d'un rectangle, est découpé en sous- 
rectangles, suivant une suite géométrique décroissante (Fig. 1.16). Chaque rectangle 
contenant au moins un point de l'objet est décompté. 

La mesure doit néanmoins rester représentative sur l'ensemble de l'espace d'étude. 
Cette condition peut difficilement être satisfaite pour certains objets (les zones 
d'échantillonnage ne sont pas toutes rectangulaires!). Aussi le calcul de la dimension 
fractale est-il biaisé pour des rectangles dont la surface dépasse le domaine de validité de 
la mesure. Une solution consiste à limiter l'espace d'étude à l'intérieur de l'objet. Mais 
cette procédure peut réduire sensiblement le nombre d'échantillons et diminuer en 
conséquence la précision du calcul. 

5.2.2.2. Longueur de la côte de la Bretagne (méthode des compas, sphère unité de 
dimension 1 ) 

On peut illustrer la méthode des compas avec l'exemple de la mesure d'une courbe très 
irrégulière, qui peut dans une certaine plage d'échelle, être assimilée au problème très 
concret de la mesure d'une côte, d'une frontière ou d'un périmètre (Fig. 1.17). 

En appliquant (1) à une longueur étalon e\, il vient: 

M E (£i) = Ni£i 

Si on réitère la mesure avec une longueur étalon Si inférieure à E\, nous obtenons: 
Me (£2) = N 2 £ 2 avec M E (£2) > M Ê fo) 
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La longueur mesurée dépend donc de l'unité de mesure choisie et croit indéfiniment 
lorsque e-K). En 1961, Richardson étudia les variations de longueur de diverses côtes et 
constata que la longueur variait suivant une loi de puissance de type: 

M E (e) = N(e)e~e-P 
Il vient d'après (2) que la dimension fractale des côtes naturelles est égale à (1+ p). 

52.2.3. Mesure de la variation de masse 

La masse d'un objet de densité constante (la masse désigne ici sous un terme très 
général le nombre de points appartenant à l'objet) centré en un point P et de rayon R, 
varie comme Rd, si l'objet est de dimension d. Les propriétés d'auto- similarité ne 
permettent pas aux objets fractals de suivre cette règle. Toutefois, on peut montrer que la 
masse de ces derniers varie comme: 

M(PJl)~R Dm 

où £> m , définissant la dimension de masse de l'objet, est équivalente à la dimension de 
Hausdorff-Besicovitch. On peut également écrire: 

M(PJi)~e d (R/ép m 

e représentant l'échelle de résolution. La dimension fractale s'obtient dès lors en suivant 
l'évolution de M(PJi) en fonction de R, par une méthode analogue à celle des 
recouvrements. Pour obtenir une mesure précise il est nécessaire d'effectuer une 
moyenne de D m sur un grand nombre de points P. Ces derniers doivent être choisis 
judicieusement afin d'éviter les effets de bord (Fig. 1.18). 

5.2.3. Commentaires 

Dans le cas de fractals théoriques, la valeur de D (pente dans une représentation 
bilogarithmique) est constante quelle que soit l'échelle. Les objets naturels terrestres ne 
peuvent bien évidemment pas vérifier cette propriété. Ds sont en effet, d'une part décrits à 
une échelle globale fixe (continent, arbre, nuage, flocon), et d'autre part dépendants des 
phénomènes physiques qui les engendrent, ces derniers étant dans le meilleur des cas 
définis à l'échelle moléculaire. Les objets réels dits fractals ne vérifient donc les 
propriétés propres aux fractales théoriques que dans une plage d'échelle finie. Le 
qualificatif fractal est accordé suivant l'étendue de cette plage (minimum 2 ou 3 ordres de 
grandeur d'échelle) et par là-même représente une appréciation assez subjective. Si 
l'explication physique de la valeur de la dimension fractale reste dans bien des cas 
obscure, la détermination de seuils d'échelle est en soi une information très intéressante 
sur les changements des processus physiques qui contrôlent le système. La recherche des 
seuils de coupure justifie à elle seule, dans bien des cas, le calcul des dimensions fractales 
des objets naturels. 
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6. Statistique et lois d'échelle sur le réseau 

Ce paragraphe passe en revue l'ensemble des résultats obtenus dans le cadre de 
l'analyse statistique effectuée sur le Blavet et la Vilaine. Pour les raisons évoquées au 
paragraphe "limites de l'étude", nous nous intéresserons essentiellement aux distributions 
de longueurs et d'aires drainées, et aux relations entre ces deux variables. 

Sous la forme d'un article, le chapitre qui suit reprend en première partie les 
techniques et les notions développées précédemment. L'auteur invite le lecteur peu 
intéressé par une présentation anglophone du contexte de l'étude, à débuter la lecture au 
paragraphe "Link length distribution". 

7. Conclusion de l'analyse statistique de la Vilaine et du Blavet 

L'analyse des distributions statistiques des longueurs et des aires drainées du Blavet et 
de la Vilaine a permis d'établir un schéma d'organisation des réseaux hydrographiques 
naturels. Les réseaux étudiés présentent en effet des densités de distribution de longueurs 
et d'aires drainées similaires à celles obtenues pour des systèmes de particules gouvernés 
par un processus d'agrégation aléatoire. Les aires drainées et les longueurs de cours 
d'eau ne semblent donc pas suivre de règles particulières de connexion. La densité de 
longueur totale de cours d'eau par unité d'aire, constante sur 3 ordres de grandeur 
d'échelle, révèle une structure remarquablement homogène et confirme l'organisation 
aléatoire des systèmes hydrographiques. 

Cependant, les caractéristiques géométriques des réseaux pour des échelles de 
longueur inférieures à 1 km (cas de la Vilaine et du Blavet) ne suivent pas les statistiques 
propres aux systèmes aléatoires. Nous pouvons ainsi distinguer deux domaines d'échelle 
où les réseaux ne suivent pas les mêmes règles d'organisation. La limite entre ces deux 
domaines est marquée par une longueur critique lo bien définie. L'interprétation de lo en 
terme de processus physique ne peut pas être menée avec les seules informations dont 
nous disposons, à savoir le réseau hydrographique au 1/100 000 et le MNT au pas de 
250 m. L'analyse qualitative du relief d'un sous bassin versant du Blavet à l'aide d'un 
MNT au pas de 30 m, montre une longueur d'interfluve du même ordre de grandeur que 
/#.. La longueur caractéristique d'organisation aléatoire des réseaux hydrographiques 
pourrait être en relation avec les processus d'érosion, qui façonnent le relief sur des 
échelles de temps long. L'étude et la modélisation des processus contrôlant la croissance 
et la géométrie des réseaux hydrographiques permettra d'éclaircir ou du moins 
d'approfondir l'origine d'une échelle critique dans l'organisation de tels systèmes. 
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Scaling Relationships and Channel Networks from Two Large- 
Magnitude Watersheds in French Brittany 

A. Crave and P. Davy 

Géosciences Rennes, UPR 4661 du CNRS, Campus de Beaulieu, 35042 Rennes Cedex, France 

Abstract. We présent statistical analysis on two watersheds in French Brittanny whose 
drainage areas are about 10,000 and 2,000 km 2 . The channel system was analyzed from the 
digitized blue lines of the 1:100,000 map and from a 250-m DEM, so that the statistics is 
obtained on a population of 4,838 and 1,136 links. Link lengths follow exponential 
distribution, confirming the Markovian model of channel branching proposed by Smart 
[1968]. The departure from the exponential distribution for small lengths, that was 
extensively discussed before, results from a statistical effect due to the finite number of 
channels and forks. The Strahler topology applied on channels defines a self-similar 
organization whose similarity dimension is about 1.7, that is clearly smaller than the value 
of 2 expected for a random organization. The similarity dimension is consistent with an 
indépendant measurement of the Horton ratios of stream numbers and lengths. The 
variables defined by an upstream intégral (drainage area, mainstream length, upstream 
length) follow power-law distributions limited at large scales by a finite size effect , due to 
finite area of watersheds. A spécial emphasize was given to the exponent of the drainage 
area, a A , that was previously discussed from différent aggregating models relevant to 
channel network growth. We show that a A is consistent with 4/3, a value that was obtained 
and analytically demonstrated from directed random walk aggregating models, inspired 
from the model of Scheidegger [1967]. The drainage density and mainstream length présent 
no simple scaling with area, except at large areas where they tend to trivial values : constant 
density and square root of drainage area, respectively. Thèse asymptotic limits necessarily 
imply that the space dimension of channel networks is 2, equal to the embedding space. The 
limits are reached for drainage areas larger than 100 km 2 . For smaller areas, the asymptotic 
limit represents either a lower bound (drainage density) or an upper bound (mainstream 
length) of the distributions. Because the fluctuations of the drainage density slowly 
converges to a finite limit, the system could be adequately described as a fat fractal, where 
the average drainage density is the sum of a constant plus a fluctuation decreasing as a 
power law with integrating area. A fat fractal hypothesis could explain why the similarity 
dimension is not equal to the fractal capacity dimension, as it is for thin fractals. The 
physical conséquences are not still really understood, but we draw an analogy with directed 
aggregating system where the growth process involves both stochastic and deterministic 
growth. Thèse models are known to be fat fractals, and the deterministic growth, which 
constitutes a fondamental ingrédient of thèse models, could be attributed in river Systems to 
the rôle of terrestrial gravi ty. 
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Introduction 

The quantification of channel networks geometry has for main goals to emphasize the underlying 
organization of this apparently complex System, and to define physical processes controlling the 
évolution of channel and topography. The relationships between physical processes and geometry of 
natural structures are not straightforward; it basically requires elementary tests such as random or 
non-random organization and the range of scale for which such organizations are valid. An indirect 
but convenient way is to validate models using empirical laws that are known to describe natural 
Systems. Since the works of Horton [1945], many expérimental measures have shown that channel 
geometries follow empirical laws. Horton's laws, Hack relation [Hack, 1957], Melton law [Melton, 
1958], fractal dimension of channel networks [Mandelbrot, 1983; Tarboton et al., 1988, 1990; La 
Barbera and Rosso, 1989, 1990] are examples of laws that are commonly used to validate models 
of channel-network growth and érosion. Most of thèse laws describe the relationships between the 
différent scales of a channel System, and seem to dérive from very basic constraints applying on 
channel networks. However they are poor indicator of the détails of the processes involved 
[Hemlinger et al, 1993; Kirchner, 1993, 1994a, b; Masek and Turcotte, 1993; Troutman and 
Karlinger, 1994]. Besides, thèse laws are predicted by most models, even built on différent 
concepts and theroritical basements (see, for example, the random walk theory [Leopold and 
Langbein, 1962; Meakin et al., 1991], or models using various deterministic rules [Roth et al, 
1989; Howard, 1990; Willgoose et al., 1990, 1991; Rodriguez Iturbe et al., 1992a, b; Rinaldo et al, 
1992; Rigon et al, 1993]). Some laws (e.g. some Horton's laws, the upstream length distribution) 
seems to dérive from very gênerai aggregation constraints applying on any channel network 
[Takayasu et ai, 1988; Rodriguez Iturbe et al., 1992a], others (e.g., the fractal dimension) are 
gênerai indicators of the interactions between instabitilites during the phase of growth. 

Before going further into models, it is fundamental to state for the validity of the différent laws 
on real channel networks. Moreover the geometrical analysis of natural Systems has to go back to 
completementary basic but detailed geometrical analysis, based on statistically relevant set of data. 
This was already mentioned by Shreve [1969] who wrote that "clearly field measurements in real 
channel networks are badly needed". 

To perform reliable statistics of channel geometry, we have chosen natural objects that fulfill 
three main conditions : 1) a significant population of links, 2) the consistency of the hydraulic 
System, and 3) a wide range of reliable informations. Point 2 implies to treat indépendant 
watersheds separately. This is important as long as the variability between différent Systems is not 
fixed; for instance, one could expect the scale of the largest watershed to fix some geometrical 
properties. Our choice was two indépendant watersheds of large-magnitude, located in the French 
Brittany, and whose outlets are in the Atlantic Océan: the Blavet (2,060 km 2 , magnitude=568), and 
the Vilaine (10,500 km 2 , magnitude=2419). Statistics are presented for each watershed separately. 
The comparison between both networks helps us to fix the effects of the watershed dimension. 
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Géométrie features as length and drainage-area distributions are analyzed on a large population, 
greater than in previous works [Smart, 1968; Krumbein and Shreve, 1970; Abrahams and Miller, 
1982]. We especially focus on the research of the best fitted law, and on the détermination of the 
cutoff scales. 

Data set and method 

Data set and geographical setting 

Networks studied in this work belong to the Vilaine basin (10,500 km 2 ) and the Blavet basin 
(2,060 km 2 ) which are both located in French Brittany (Figure 1). The basement is mainly 
constituted by schists (Vilaine) or mixed granités and schists (Blavet). To the South of the 
Armorican Massif, the topography is partly controlled by a séries of Hercynian shear zones, 
oriented at NI 00°. Thèse accidents présent a slight récent tectonic activity, marked by moderate 
earthquakes whose magnitudes never exceed 4.5 on the Richter scale. The précipitation on both 
watersheds is about 600-700 mm.a" 1 . 

For both watersheds, we used two types of data: 1) the Vilaine and Blavet digitized networks, 2) 
the D.E.M. of the French Brittany giving élévations on a 250-m square grid. The source of the 
digitized networks is the blue Unes on the 1:100,000 "Institut Géographique National" maps, 
resolution that fixes the définition of headstream position. We are able to check géométrie features 
as length of channel links, on several orders of scale, with an accuracy of about 100 m. We have 
chosen the 1:100,000 map in relation to the scale of the studied area, and as a compromise between 
technical problems (in particular, the digitalization) and accuracy. Apart from the inhérent accuracy 
on channel location that is estimated at about 100 m, the most important problem is the définition 
of the smallest détectable features and of their real order. A comparison with the 1:25,000 map and 
field observations shows that the most of the small channels are detected with the 1:100, 000 map 
and that only the length of the order- 1 channels is affected by the précision of the map. However, 
the précise location of headstream is a quite unsolvable problem. Headstreams présent in fact strong 
variability of extension with climatic conditions and can be hardly defined even in fields, when 
sources are diffused saturated areas or graduated channels [Montgomery and Dietrich, 1988, 1989; 
Kirkby, 1988]. We have verified on some local points that the 1:100, 000 map gives a rather good 
average of permanent headstream location. 

Numerical processing 

The channel network topology, is determined with the assumption that channel networks are 
perfect tree networks, without divergences in downstream direction. This implies that river path is 
unique, between any point of the channel network and the outlet. In most studies, the "unique path 
assumption" is widely used as an implicit assumption of channel networks and the topological laws, 
like Horton's laws, are only defined under this condition. Natural channel networks, like the Vilaine 
and Blavet, show some proportion of downstream divergences, mostly located in very flat areas. 
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But the low proportion of such divergences (less than 2 %) within the Blavet and Vilaine structures, 
allows us to use the unique path assumption for both networks. In that cases, the existent 
downstream dérivations was removed before any numeric calculation by a "cleaning" numerical 
routine that eut one of the two tributaries of each downstream dérivation and transforms the initial 
digitized networks into "perfect" trees. It is then possible to build the "family tree" of ail links and 
to recombine them following topologic rules. An other routine calculâtes for each link its length 
and its total upstream length. Links are recorded in a vector lines form, so length and upstream 
length can be calculate without artifacts inhérent to bitmaps, and due to the pixel size. 

Drainage area can be calculated if flow directions are known at any DEM's point. To satisfy the 
unique downstream path condition, only one flow direction was defined. The local steepest-slope 
criteria is a practical mean to détermine a flow direction [Mark, 1983; Mark et al, 1984; Band, 
1986, 1989]. It is implicitly assumed that flow processes are exclusively controlled by slope 
parameters. Such assumption is realistic for land-systems where surface runoff is the dominant flow 
process. But in many cases, especially at micro-scale, water transfer on hillslope cannot be 
simplified to this simple model [Dunne, 1980; Kirkby, 1988]. In French Brittany, most of the water 
flow is under subsurface form, expect during sparse intensive storm events when some runoff 
processes are observed. The slope criteria only applies if the bedrock surface is parallel to the 
topographie surface. However, this assumption seems reasonable, considering the 250 m resolution 
of our DEM. Hillslopes are, for most of them, cover with only 3 or 4 pixels. Meanwhile the 
estimation of drainage area for hillslope is inaccurate. 

Gridded élévation data usually présent some lower points where ail neighbouring grid points are 
more elevated. Even if natural network could display such configurations, it is most of the time due 
to the artifacts of the interpolation procédure. Published numerical procédures adjust the drainage 
direction for this configuration cases, and ensure the continuity of the drainage network [Band, 
1986, 1989; Morris and Heerdegen, 1988]. We also used information from the digitized channel 
networks to improve the local drainage direction. Thèse computations are carried out by the 
MNTSurf software (unpublished), from the "L'atelier de spatialisation numérique" from Ecole 
Nationale Supérieure Agronomique de Rennes. MNTSurf has been modified in a way to take into 
account digitized networks. 

Statistical properties of channel networks 

Ail the entities and terminology reported in this section are defined in the works of Shreve [1966] 
and Abrahams [1984]. "Link" refers to a channel between two forks, "stream" to a channel entity 
defined by the topological ordination from Strahler scheme. For each distribution length we 
presented: 1) the density distribution n(l), where n(l)dl is the number of links of length between / 
and l+dl, and 2) the cumulative distribution C(/) the number of link of length larger than /, i. e.: 




(1) 



L 
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Statistical studies on populations usually calculate C(Z) to underline the main tendency of the 
statistical set. To improve the analysis we also calculate n(Z) for two reasons: 1) fluctuations of «(/) 
give the order of précision of statistical parameters, and 2) n(l) are less sensitive to finite-size 
effects at large lengths than C(Z). A separate calculation, and a comparative study of «(/) and C(Z) is 
a great helpfull method to define the more pertinent theoretical fitting distribution. Indeed the ratio 
between the density distribution and the cumulative distribution n(l)/C(l) takes characteristic forms 
for typical functions. It should be for example independent of length for exponential function and 
should vary as the inverse of the length, l~ l , for power laws. The method to calculate density 
distribution is original and developed in Davy [1993]. This procédure calculâtes, by exploring ail 
the range of link length, the length interval dZ for which the distribution n(l)=N/dl présents the 
lowest fluctuation, with N the number of links in each interval. To account for the large variations 
of the length scales (from 0.1 to 10,000 km), n{ï) and C(Z) are calculated for both a linear variation 
of lengths (constant différences between two consécutive lengths) and a logarithmic variation 
(constant ratio between two consécutive lengths). 

Once a distribution model is established by the comparison between «(/) and C(l), we calculate 
the model parameters using the smoothed distribution C(Z), and we estimate parameter error bars 
using n{l). 

Link-length distribution 

Since the beginning of studies on channel geometry, nature of link length distribution is debated 
[Abrahams, 1984]. The search of the best-fitting law is as important as the collection of data since it 
is the support of any hypothesis on processes. But observations on several channel networks did not 
emphasize an unique distribution law [Smart, 1972; Dunkerley, 1977; Abrahams and Miller, 1982] 
since lognormal, gamma, mixed-gamma functions have been used as the best fits of the link length 
distribution. Also exterior and interior links were considered to follow différent distributions, but 
the relationships between thèse two populations of links remain an open question [Schumm, 1956; 
Jarvis, 1976]. 

Density distribution of link length is reported for both Blavet and Vilaine networks in Figure 2. 
To compare the results from the both networks, n(l) is normalized by the total number of links for 
each network. It shows that link length distribution follows the same trend whatever the network. 
n(l) présent two domains of variation with an increase from 0 to 1 km, and an exponentially 
decrease for larger length. For length greater than 1 km, n(Z) varies as erWo, where Z 0 is the 
characteristic length of the distribution, equal to 1±0.1 km. Neither lognormal function nor power 
law are appropriate to fit the data. This could be emphasize by the plot of the ratio n(Z)/C(Z) which 
remains constant for length larger than 1 km. Note that interior and exterior link length distributions 
follow exactly the same trend (Figure 3). This point was debated by Abrahams [1972] and 
Abrahams and Miller [1982] on smaller watersheds. This study show that there is no statistical 
reason to make a distinction between both types of links. 
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The question of small length distribution is of great importance since it conditions the search of a 
best-fitting model. Différent interprétations can be proposed, depending on the necessity to model 
small length distribution that is here below 1 km. If small lengths really belong to the distribution, 
we have to look for a distribution that ensure the existence of a peak at finite length (for instance, 
lognormal or mixed-gamma function discussed in Abrahams and Miller [1982]). If there are 
physical or observation-accuracy reasons for removing small links from the distribution, the 
exponential function is a very efficient fit which has a small length cut-off at about 1 km. The need 
of interpretating small length distribution is a récurrent problem for natural data; let us quote 
Einstein and Baecher [1983] who pointed out that lognormal fits for fault length distributions only 
quantify the lack of observation of small structures ! 

To interpret small length distribution, let' s first notice that there exists a limited number of links, 
and subsequently of forks. This must introduce a limit of the link distribution, whose simplest 
model would be the distances between nearest fork neighbours. If we assume an average density of 
forks, ^Fork=0.7±0.2 km" 2 (this point will be addressed in a next section), this simple model leads to 
a Rayleigh distribution that is: 

f 

n(l)=2 iV ap p k d app l exp(-7C d app l 2 ), (2) 

where iV ap p and d app are the number and density of links that follows this distribution. If we assume 
that a portion a of the total number of links follows the nearest-neighbourgh distribution (N app =ccN 
and J ap p=c^Fork), we find the distribution of small length is well fitted for a between 0.5 and 0.7 
(Figure 2). The Rayleigh function is defined by one main parameter which is équivalent to a length 
scale : 



This length Z R is about 1 km, that is identical to the length scale l 0 derived from the exponential 
function at large lengths. It is not surprising since both distributions are basically governed by the 
density of links. 

The nearest-neighbourg model however générâtes loops between forks, and is thus inconsistent 
with a tree organization. A simulation of a tree-like organization based on a random distribution of 
forks shows that, for lengths greater than the peak of the Rayleigh distribution (i.e. the maximum 
probability), the distribution becomes an exponential function. 

Our results are consistent with the "Random Topological Model", proposed by Shreve [1966] 
who postulated that channels are topologically random. It means that the probability of finding a 
junction is equivalently distributed within the network (i.e. a Markovian process of branching), as 
expected by Smart [1968]. Since this work, no set of data on expérimental channel networks 
supports clearly thèse theoretical results [Abrahams, 1984]. The main problem that we raise is the 
pertinence of small structures that are generally considered belonging to the distribution model of 
link lengths. The departure from the exponential distribution for small links is intrinsically related 
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(3) 
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to the fini te number of links that introduces a lower-cutoff, below which data are not représentative 
to the distribution. The lower cut-off and the large-length distribution are both an unique 
characteristic length scale, / 0 ~1 km, similar for both watersheds. It is clear that l 0 is defined by the 
density of links, which is related to our capability to describe small rivers, i.e. by the resolution of 
the blue lines. This length in Brittanny could represent more a topographie length scale than an 
map-operator scale. An example from an order-1 subcatchement (Figure 4, obtained from a DEM 
with a resolution of 27 m) show that the topography roughness is obtained with wavelengths larger 
than 1 km. This length scale can be observed in many places as a significant cutoff of the 
topographie spectrum of the breton landscape. For a long time, the existence of a cutoff length has 
been proposed to explain the lack of incisions up to di vides [Horton, 1945; Montgomery and 
Dietrich, 1988, 1989], that is a real conséquence of the erosive process. However we do not want 
to stress too much this observation that requires a more gênerai analysis of the topographie 
spectrum. 

Length distribution of Stahler streams 

To better describe the tree organization, we use the Strahler ordering topology that defines a 
hierarchy between streams [Strahler, 1952]. The comparison between streams of same ordering 
(quoted co) allows us to refine the description of the tree structure. The length distribution of 
streams of order 2 and 3 are shown in Figure 5. Because of the reduced number of data, the 
distributions présent large fluctuations; but an exponential fonction is always consistent with the 
data with a characteristic length scale, l m , that increases in a geometrical progression with the order 
co (Figure 6). The number of streams, N w , also decreases in a geometrical progression with Où. This 
resuit is consistent with Horton's laws [Horton, 1945], which state that the ratios Rb and R\ (Rb=N (û / 
Nfn+i and i?i=</> aM .i/<Z> (0 , with </>£,, denoting the average length of order co) are constant whatever 
the order co. Note that, if we assume that the exponential distribution of link length is représentative 
to the tree organization, R\ is better calculated by using the length scale l a . In theory, the average 
and the characteristic length scale are equal for an exponential distribution (also is the standard 
déviation). In practice, the détermination of is not sensitive to the population of non- 
representative small streams; but the average <l> a is. 

We have calculated Rb and i?i from the slope of the linear régression of log(A ? û) ) and \og(l w ) versus 
G). R\ is equal to 2.4+0.25 and 2.1±0.25 for the Vilaine and Blavet networks, respectively. Rb is 
equal to 4.2±0.1 and 4.15±0.1 for the Vilaine and Blavet networks, respectively (Table 1). Previous 
studies on stream networks have shown that Rb and R\ are restricted to narrow ranges: 3<i?b<5 and 
1.5</?i<3 with modal values equal to 4 and 2, for Rb and R\ respectively. For both studied networks, 
thèse Horton's ratios are closed to the modal values. 

Another interesting distribution is the total stream length distribution which mixes streams of ail 
orders. This distribution is différent from the link length distribution because a Strahler stream can 
be formed by several links. This distribution was used by Tarboton et al. [1988] to define the fractal 
dimension of channel networks. The stream distribution follows a power law trend t a s (Figure 7) 
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for stream lengths larger than 1 km. The exponent as is equal to 1.65±0.05 and 1.7±0.1 for the 
Vilaine and Blavet networks, respectively. The main source of uncertainties is the choice of the 
relevant length scales. The estimation is better for the larger basin because the dynamics of length is 
larger. 

Because the number and lengths of order G) are in a geometrical progression, one can show that 
the exponent as must be related to R\> and R\ by the following expression: 

a s =log(i? b )/log(fl,), (4) 

The démonstration was proposed by Tarboton et al. [1988] if streams of order co have a constant 
length. The resuit also applies to a sum of exponential distributions if their characteristic length 
scales /a, and number of streams N m are in a geometrical progression (Figure 8). The exponent as 
quantifies the hierarchical composition of the river network. In a sensé, it looks like a fractal 
dimension, called the similarity dimension by Feder et al. [1989], by analogy to deterministic self- 
similar fractals. For the diffusion-limited aggregation and viscous fingering clusters, as is 
équivalent with other fractal dimensions (capacity dimension, mass dimension, the radius-of- 
gyration dimension, ...) [Feder et al, 1989]. But, the existence of similarity dimension, based on a 
topological ordering, does not imply necessarily that the network itself is fractal. This point will be 
discussed later. 

Taking the most probable values, i?b=4 and R\=2, which correspond to a topologically random 
organization [Shreve, 1966], we should have obtained as=2. The exponent measured for the Vilaine 
and Blavet networks (1.65 and 1.7, Table 1) show that both networks are far to satisfy the random 
model. This values are consistent with the Rb and R\ ratios measured independantly (Figure 6). 
Moreover, because the accuracy on the exponent as is relatively good, it necessarly implies Ri ratios 
significantly larger than 2 for both networks (2.3-2.5 for the Vilaine network, and 2.1-2.4 for the 
Blavet network). 

A complète characterization of the System would have required to calculate the stream area ratio, 
i?A, that fixes the ratio between average areas drained by each successive order stream. This ratio 
was initially proposed by Schumm [1956], and its implication on the détermination of scale 
dimension is discussed by Rosso et al. [1991] and Beer and Borgas [1993]. However, as for other 
ratios, R A is accurately defined only for low-order drainage basins where the number of éléments is 
sufficient; and the pixel resolution of our DEM is not appropriate to perform this calculation. Ra is 
expected to be larger than Rb and equal to i?i 2 for an idéal self-similar Horton System [Beer and 
Borgas, 1993]. R A is equal to Rb if streams are space filling [Rosso étal, 1991]. 

Upstream length, drainage area and mainstream length distributions 

The upstream parameters (upstream length, drainage area or mainstream length) are 
représentative to the aggregation properties of the network. The drainage area, A, and upstream 
length, Z, are the integrate of areas and river lengths that drains into a given point of the watershed. 
The mainstream length, lu, is the longest channel length that connect a point of the network to a 
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source, and it was considered as a measure of the flow distance in watersheds. Thèse distributions 
are assumed to follow power-law relationships, regarding previous analysis on natural streams 
[Rodriguez Iturbe et al, 1992a] or aggregation models [Leopold and Langbein, 1962; Takayasu et 
al, 1988; Mandelbrot, 1989; Meakin et al, 1991; Masek and Turcotte, 1993; de Vries et ai, 1994]. 
The distributions are estimated from a set of variables (A, Z or lu) calculated at the downstream 
termination of each link. We made no distinction between exterior and interior links because both 
link length distributions are similar. 

AH thèse distributions (figures 9, 10, 11) shows a lower cut-off at about 1 km for lengths and 
1 km for areas, consistent with the finite extension of the tree network, as it is discussed above. For 
several order of scale ranges, the density distributions show remarkable power-law variations with 
exponents roughly around 1.5. Further interesting information, concerning finite size effects, can be 
obtained from the combined use of the density and cumulative distribution. Actually, the 
cumulative distribution is clearly no more a power law at large lengths or areas while the density 
distribution is. If the density distribution is really a power law, this différence must be related to the 
upper bound, s max , of the intégral : 



T f s 

>= )n\ 



C(s) = J n(s) ds 



V^max J 



-1, (5) 



where s is the scale variable of the distribution (area or length), and s max a finite size effect related 
to the extension of the watershed. This large-size effect was derived mathematically from a Horton 
system by La Barbera andRoth [1994]. 

The cumulative distribution départs from the power law, even at small scales if the exponant a of 
the distribution is not much larger than the critical value of 1, a value for which C(s) is dépendant 
on s m ax at ail scales (C(s)~\n(s/s max ) for a=\). This finite-size effect can induce large overestimation 
of the exponent a, if it is directly measured from the logarithmic slope of the cumulative 
distribution. To illustrate this point, let's calculate the exponent a from the logarithmic slope of 
C(s): 

dln(C) \-a 



dlns " f s 
1- 

V^max J 



(6) 



If a précision £ is required for the expOnent a, only scales smaller than s max *£ 1/< ' a ' l) must be 
considered. For instance, a précision of about 20% (£=0.2) is obtained for scales smaller than 
Smax/25 if a=\.5, or 5 max /215 if a=\3. Considering the range of relevant scales for our data, this 
constraint gets difficult a reliable détermination of a from the logarithmic slope of C(s), except if 
s m ax is explicitely incorporated in the régression fit. This example also illustrâtes the potential 
efficiency of a combined use of the density and cumulative distribution. 

The upstream length distribution is only weakly sensitive to finite size effects (Figure 9). For the 
Blavet network besides, both density and cumulative distributions follow power-law trends at ail 
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scales. The measured exponents, az, of the upstream length distribution are 1.42 and 1.45 for the 
Vilaine and Blavet networks, respectively (Table 1). The error bar on thèse exponents is about 0.02, 
so that thèse two values are not significantly différent. 

The drainage area distribution shows finite-size effects at about 6500 km 2 and 2400 km 2 for the 
Vilaine and Blavet basins respectively (Figure 10). The estimate of thèse values is rather poor since 
it is mostly constrained by the end of the distribution where large statistical fluctuations are présent. 
Considering error bars on A max , the finite-size effects is well interpreted as the whole surface of 
both watersheds. The exponents of the drainage area distribution, a A , are 1.3 and 1.35 for the 
Vilaine and Blavet basins, respectively (Table 1). Once more, thèse values are not significantly 
différent, considering the error bar of about 0.05. Rodriguez Iturbe et al. [1992a] have measured 
rather larger exponents for five basins from North America, with average values about 1.4 (0.4 for 
the cumulative distribution). We could have derived similar values from our data, if large-size 
effects are not considered (see équation 6). 

The mainstream length distribution shows very important finite size effects of about 160 and 
140 km for the Vilaine and Blavet networks, respectively (Figure 11). The exponent of the 
distribution, a M , is well defined for the Vilaine basin (<z M =l-43±0.02). For the Blavet network, the 
régression give a larger exponent of about 1.7. However, the fit is not accurate for lengths between 
1 and 5 km. A smaller value (1.5, dashed line in Figure 11) could be obtained if lengths between 1 
and 5 km are considered statistically significant. 

The interprétation of thèse exponents in terms of channel characteristics was discussed through 
growth models of channels [Meakin et al, 1991] or aggregation properties of simple trees 
[Takayasu et al, 1988]. The simplest meaningful values for the exponent a A (which is the only 
exponent extensively discussed) are 3/2 for compact river basins filling a planar région, and whose 
width and length scale similarly. Another trivial value of 4/3 was proposed for self-affine drainage 
basins whose width scales as the square root of the height. Such scaling is intrinsically obtained for 
models where the displacement is random in one direction and constant in the other, as the random 
walk's model of Scheidegger [1967], the. random aggregation model of Takayasu et al [1988] and 
the directed self-avoiding random walk of Meakin et al. [1991]. Intermediate value of 1.39 was 
obtained from a self-avoiding random walk by Meakin et al. [1991]. The measured drainage 
exponents a A are found to be intermediate between the two extrême values of 4/3 and 3/2. 
Rodriguez-Iturbe et al. [1992a] and La Barbera and Roth [1994] found exponents around 1.4; the 
former proposed that the departure from the trivial exponent 1.5 is related to the fractal dimension 
of the flow path (here the mainstream). Our data (Figure 10) are clearly inconsistent with an 
exponent of 3/2, and consistent with an exponent of 4/3. But, considering error bars, mainly due to 
the choice of the représentative area range for fitting, exponents as large as 1.4 cannot be definitely 
excluded. Note that the model of Takayasu et al. [1988] makes no différence between the drainage 
exponent a A and the upstream exponent a z , which is equal to 1.45 for both watersheds. 
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Density properties 

The question of the spatial organization of channels was adressed from the scale dependency of 
channel density that relates the total upstream channel length, Z, to the total upstream drainage area, 
A: 

D A = Z/A. (7) 

This measure is exempt from any topological considération, but obviously dépends on the 
accuracy on the description of the network. If the density D A varies with the area of the watersheds, 
the channel network is a fractal object with non-trivial spatial property. This point is rather 
important since the fractal dimension of the channel network can give important constraints on the 
underlying mechanisms of channel growth (see for instance [Stark, 1991; Meakin et al., 1991; 
Masek and Turcotte, 1993]). The previous analysis of the distributions of drainage area and 
upstream length have shown that both can be described with a power law whose exponents are 
significantly larger for the upstream length (az~1.45) than for the drainage area (aA~1.35). Thus we 
expect a non-linear relationship between Z and A (Z-A&i), with an exponent a,j about 0.7 to 0.9, 
following : 

l-a A 
l-a z 

From a technical point of view, an estimate of the drainage density requires an adjustement 
between the digitized network and the DEM. This adjustement can give abnormal parameters for 
some points of the network, because 1) the network is defined at ±1 pixel, and 2) the pixel size of 
the DEM can lead to abnormal drainage configurations when compared to the real case deduced 
from the digitized network. When first-order catchment are removed, it remains about 10% of 
abnormal values. Most of them were easily détectable and removed, but about 1% of the data set 
presented in Figure 12 still présent abnormal values, certainly due to this technical problem. 

The relationship between drainage density and drainage area is presented in Figure 12. The 
average density, marked in the figure by the largest concentration of points (bold arrow), is clearly 
decreasing for areas less than 100 km 2 . For larger areas, the density is almost constant around a 
value of about 0.7 km" 1 for both Vilaine and Blavet basins. If a numerical fit with a power law is 
performed on the total set of data, it gives an exponent of about -0.2 to -0.15 that is not really 
représentative to the trend of the average density which rather shows a decrease of the average 
drainage density with area. Similar trends were observed by La Barbera and Rosso [1990] for 60 
subbasins of Alto Liri basins in Southern Italy. Thus clearly the drainage density varies with scale, 
at least up to 100 km 2 , but without fractal properties. 

Because of technical problem discussed above, a précise analysis of the data shown in Figure 12 
could be disputable. It is possible to define quite precisely the évolution of the drainage density at 
différent scales without constraints of basin boundaries, i.e. in an arbitrary area of constant form. 
This method is similar to the Mass method classicaly used to calculate the cluster dimension or 
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mass dimension of fractal objects [Feder, 1988]. The drainage density d(Â) and its standard 
déviation Oa(A) was calculated for 100 points, chosen at random on the network, and far enough 
from the boundary of the watershed. The average mass dimension (open circles in figure 13) shows 
the same trend than in Figure 12. The only différence stays at very large areas (>2 10 km ) when 
the measurement surface is partly outside the limits of the basin. Thèse large areas are thus not 
considered in the following discussion. Also small areas (<1 km ) are not représentative since they 
only include one link on average. 

In the représentative area range, the lower bound of the drainage density distribution (d-cf) shows 
a remarkable constant value of about 0.7 km" 1 , that is consistent with the asymptotic density 
deduced from Figure 12. The standard déviation Od(A) shows two remarkable trends in this area 
range. For small areas (1<A<100 km 2 ), it decreases as a power law with an exponent of about 
«d=-0.37±0.02. For areas larger than 100 km 2 , cis constant around 0.1 km" 1 . The évolution of the 
average density with area in the range 1 to 100 km should be modelled as : 

</=<L+y d A-Od, (8) 

where and % are two constant values, représentative of thè basin. This type of scaling is 
described for objects so called fat fractals. They differ from a pure fractal object by a non-nil 
density for large Systems, A— >°<> [Liang and Kadanoff, 1985]. 

Similar trends were observed on the variation of the fork density with area (Figure 14). The 
average density of fork d F and its standard deviaton <x F , significantly decreases for areas smaller 
than 100 km 2 . The standard déviation Of exhibits the same power-law trend than for the drainage 
density, with still an exponent of about -0.35. The lower bound of the fork density is constant over 
3 orders of area magnitude. 

This analysis is also appropriate to the relationship between the mainstream length Zm and the 
drainage area A, introduced by Hack [1957], which was discussed in many papers as a main self- 
similar property of basin subcatchment [Hack, 1957; Gray, 1961; Mandelbrot, 1983; Mesa and 
Gupta, 1987; Hjelmfelt, 1988; Tarboton et al, 1988; Robert and Roy, 1990; Rosso et al, 1991]. 
Following Hack, lu and A should be related by a power law such as : 

/m~A«h (9), 

where an generally départs from the trivial value of 0.5. A direct measurement of Cfa from a linear 
régression of log(/M) versus log(A) leads to a value of about 0.65 (Figure 15), that is merely 

2 2 

représentative of the scattered data at intermediate areas (<100km ). For large areas (>100km ), 
the data are best-fitted by an exponent On of 0.5, which represents a trivial scaling of the main 
channel (z M ~VÂ). This trivial asymptotic limit was already recognized from large-basin 
measurements [Mueller, 1973], or models of Markovian branching [Mesa and Gupta, 1987]. Below 
this critical area, this large-scale fit constitutes an upper bound of the fluctuations observed on the 
data plot (Figure 15). We thus wonder if the Hack' s relationship would be better described by an 
équation analogous to (8) : 
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i M ~4Â-rA p . (io) 

A complète characterization of équation (10) would require a more précise analysis to be 
completely validated. 

The scale dependency on both density and mainstream length strongly suggests that the channel 
network System cannot be described by a simple fractal geometry. The main argument is that, for 
both variables, it exists an asymptotic trivial limit for the dimension of the System. At intermediate 
areas, data show strong fluctuations that were barely incorporated in previous descriptions. We first 
state that thèse fluctuations can contain an interesting information about the channel organization 
processes. If our analysis is correct, the scale réduction of fluctuations can be described by a power 
law for areas between 1 and 100 km . The other interesting property of thèse Systems is that the 
trivial asymptotic limit represents a lower bound (for drainage density) or an upper bound (for 
mainstream length) of the whole distribution. This type of behaviour, where the average value is the 
sum of a constant term and a power law, was mentionned by Mandelbrot [1983] as non-lacunar 
fractals whose dimension is equal to the embedding dimension, and by F armer [1984] as "fat 
fractals". This type of behaviour was suggested for the scaling of the drainage basin (total area 
minus the surface area of river channels) by Karlinger and Troutman [1992]. 

More interesting for our analysis is the récognition of fat-fractal properties for a ballistic 
aggregation model (or directed walk model) [Liang and Kadanoff, 1985; Meakin, 1991; Vicsek, 
1989], where the particles move in a random way in one direction (say x) and with a constant 
velocity in the other direction (say y), before sticking to a growing cluster. Such growth model 
produces a fan-like structure (see figure 7.4 p. 195 in [Vicsek, 1989]), whose average density follows 
équation (8). Because of the displacement condition on y, the fan is a regular cone, meaning that its 
average width is proportional to the structure radius. This regular geometry forces the System 
dimension to be the embedding space dimension (this was demonstrated in Vicsek [1989] by using 
an argument on the growth rate of the radius). However the system is locally inhomogeneous, and 
the average density converges slowly with scale to a constant value, such as d-d^-f^ [Liang and 
Kadanoff, 1985]. The exponent f5, that is twice our exponent 0Cd, was obtained from simple 
aggregation models with non-universal values (i.e. depending on the lattice topology, cone angle, 
etc..) varying between 0.55 and 0.9 [Liang and Kadanoff, 1985], that is consistent with the 
exponent measured from channel network (/?=2ad~ 0.74). As for the directed-walk models, channel 
growth cannot be considered as a pure aleatory aggregating System because the average flow is 
driven by the terrestrial gravity along the average basin slope. The new results that we have derived 
from the observation of natural system should lead to push further the comparison between the 
différent channel growth model, in the same spirit of the work presented by Meakin et al. [1991]. 

Discussion about the Nature of Channel Network. 

The difficulty to interpret geometrical properties of channel network is to connect two apparently 
différent problems that are the topology of the branching, and the spatial distribution of the 
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network. At first, we présent separately the results which concern the "tree" organization of the 
channel network, and those which concern the spatial organization of the channel network. 
The main results for the tree organization, i.e. obtained in the referential of the tree, are : 

1. Indvidual links or streams follow exponential distribution. The departures from the 
exponential distribution is only a statistical effect due to the finite number of links. 

2. As a conséquence, the branching is almost equally probable everywhere on the tree, as 
expected by Smart [1968]. 

3. The topological ordering of Strahler defines a self-similar System, whose similarity dimension 
is clearly less than 2 (about 1.6-1.7). The similarity dimension is completely defined by the 
Horton ratios, but the most accurate characterisation is obtained from the distribution of ail 
Strahler stream lengths. 

4. Variables defined by an upstream intégral (drainage area, mainstream length, upstream length) 
follows power-law distribution whose exponents contain informations about the aggregating 
process (Table 1). 

The results obtained in the plane referential, and concerning the spatial organization of channels, 
are : 

1. Both density and mainstream length have a trivial asymptotic limit corresponding to a 
constant density and square root of areas, respectively. This necessarily implies a space 
dimension of 2 for the channel network, equal to the embedding space. The asymptotic limit is 
reached for drainage areas larger than 100 km 2 . For smaller areas, the asymptotic limit 
represents either a lower bound (drainage density) or an upper bound (mainstream length) of 
the distributions. 

2. Because the fluctuations of the drainage density slowly converges to the finite limit, the 
System could be adequately described as a fat fractal, a model which predict the density 
fluctuations to decrease as a power law with integrating area. 

In simple Systems, as thin deterministic fractals, there exists a direct relationship between spatial 
properties (such as a fractal dimension), and hierarchy of elemetary structures (number and length 
of streams) that composes the total network. Starting from the combining ideas that channel 
networks are Hortonian (i.e. there exists a hierarchy of entities, ordered whose length and number 
are in a geometrical progression with the Strahler ordering parameter), previous studies have related 
the parameters of this topology, mainly the Horton ratios, to spatial dimensions of the network. In 
other words, the similarity dimension which defined the hierarchy organization of the tree is 
considered as equal to the fractal (capacity, mass) dimension. This was demonstrated for DLA and 
viscous fingering structures by Feder et al. [1989]. In channel networks, the similarity dimension is 
clearly smaller than the capacity dimension. This discrepency was interpreted as a conséquence of 
the fractal nature of individual streams [Tarboton et al, 1988, 1990], or as the way the stream 
length or the order of the watershed scales with the resolution of the map [La Barbera and Rosso, 
1990; Beer and Borgas, 1993]. A more basic reason for the différence between similarity and 
capacity dimensions could be that the implicit assumption that channel network are thin fractals is 
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wrong. For both Blavet and Vilaine networks, the drainage density is better described by a fat 
fractal for areas smaller than 100 km 2 , meaning that the density has a finite limit for large areas. An 
interesting conséquence of such assessment is that informations about the spatial organization, and 
about the underlying aggregation process, could be contained in the scale dependency of 
fluctuations, rather than of average density. The critical area of 100 km 2 could represent a kind of 
limit for the aggregation process, for which the System is completely homogeneized, i.e. the natural 
random fluctuations become larger than "aggregation" fluctuations. 

The physical meaning of thèse observations is not straight, because 1) scale exponents and 
growth process are not linked by an univocal relationship, and 2) the littérature on the theoritical 
background of fat fractals is, to our knowledge, still scarce [see however Eykholt and Umberger, 
1988, 1992]. Several studies have tried to characterize the channel formation process, either from 
dynamical érosion models, or from simple aggregation models; the former models give mainly 
informations about the physics of érosion, channel geometries being used as a validation parameters 
for models, while the latter gives information about the physics of the channel organization. Ail 
thèse models base the validation procédure on the exponents a A and a M , the Horton ratios and the 
similarity dimension as of the channel network. 

The statistics on channel network, and more precisely the departure of scale exponents from 
trivial values, tend to prove some determism in the network building. The first argument is from the 
similarity dimension as, which is significantly smaller than the trivial value of 2 (equivalently the 
Horton length ratio R\ is significantly larger than 2). A second argument is from the drainage area 
exponent a A , whose departure from the trivial value of 3/2 suggests a determinism in the network 
building [de Vries et al, 1994]. The uppermost value of a A , consistent with data presented in Figure 
10, is about 1.4, with a mostly probable value of 1.35. The question now is to précise the type of 
determinism that could be consistent with our results. The very simple aggregation models, such as 
the diffusion-limited aggregation, are particularly attractive to relate a growth process to scale 
exponents. Thèse models simulate the growth of a cluster from the aggregation of walkers. 
Différent walking and sticking rules produce différent types of geometries; but there exists only a 
few universal classes, characterized by the same set of scale exponents. Because a river System is a 
kind of cluster, thèse models were often applied to understand scale relationships of river Systems 
[Leopold and Langbein, 1962; Scheidegger, 1967; Seginer, 1969; Howard, 1971; Takayasu et al, 
1988; Meakin et al., 1991; Stark, 1991; Masek and Turcotte, 1993]. Ail thèse models should be 
revisited for a full characterization of the scaling of the drainage density. 

Models based on directed random walk aggregation should lead to a fat fractal geometry, as 
predicted by Liang and Kadanoff [1985]. The simplest directed random walk [Leopold and 
Langbein, 1962; Scheidegger, 1967; Takayasu et al, 1988; Meakin et al, 1991] predict the 
drainage density distribution exponent a A to be exactly 4/3, consistent with the exponents measured 
from the Vilaine and Blavet network. The similarity dimension and Horton ratios are not 
inconsistent with the one measured from both breton catchments (i?i about 2.8 [Leopold and 
Langbein, 1962] or 2.9 [Meakin et al, 1991]), but further work is required to define more precisely 
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thèse scaling exponents. Another reason to emphasize this model is that the terrestrial gravity 
imposes a deterministic flowing condition to natural Systems, which is the necessity for rivers to 
follow the main slope of the basin. We need more arguments to validate this assessment but we 
think that models based on the combination of deterministic and random conditions could represent 
an interesting prospective way to study river network scalings. 
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Figure Captions 

Figure 1 . Geological map of the Armorican massif with the location of the Vilaine 
(A) and the Blavet (B) basins. 

Figure 2. Density distribution of link lengths for the Blavet and Vilaine networks. 
A Rayleigh function (dashed line) fits small length distribution. The Poisson 
(exponential) function fits well lengths larger than 1 km. 

Figure 3. Density distribution «(/) (black squares), cumulative distribution C(7) 
(bold lines) and exponential fit (thin solid line) for the link length distribution of 
the Vilaine network. Top, exterior links; bottom, interior links. Similar trends are 
observed for the Blavet. 

Figure 4. Topography of the Coet-Dan subcatchment (Blavet basin), obtained 
from a DEM with a resolution of 27 m. 

Figure 5. Density and cumulative distribution for stream lengths of Strahler order 
2 (left side) and 3 (right side) for the Blavet (top) and Vilaine (bottom) networks. 

Figure 6. Plot of the number of streams N a (top), and of the stream length la, 
(bottom) as a function of the Strahler order û) for the Vilaine (dark squares) and 
Blavet (open circles) networks. l m is calculated for stream orders smaller than 4, 
for statistical reasons. The fitting curves (solid line) define the Horton ratios R b 
and R\. 

Figure 7. Density (solid squares) and cumulative distribution (bold lines) of 
Strahler-stream lengths for the Vilaine (top) and Blavet (bottom) networks. The 
power-law fits are drawn with a thin solid line. 

Figure 8. Plot of Cil) (bold line) defined as the sum of exponential functions 
(dashed lines), whose characteristic length scales and multiplicative coefficients 
are in a geometrical progression, defined by ratios R\ and i?b respectively. The 
exponent derived from the slope of log(C(/)) versus log(/) is exactly equal to the 
ratio log(i? b )/log(i?i). 

Figure 9. Upstream length distribution for the Vilaine (top) and Blavet (bottom) 
networks. Same symbols as for figure 3. 

Figure 10. Upstream drainage area distribution for the Vilaine (top) and Blavet 
(bottom) networks. Same symbols as for figure 3. 

Figure 11. Mainstream length distribution for the Vilaine (top) and Blavet 
(bottom) networks. Same symbols as for figure 3. The dashed line shown for the 
Blavet network indicates a power-law fit with an exponent of -1.5 (see text). 
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Figure 12. Plot of the drainage density versus integrating drainage area for the 
Vilaine (top) and Blavet (bottom) networks. The bold arrow indicates the average 
trend. The solid line indicates the asymptotic limit. 

Figure 13. Scale dependancy of the channel density for the Vilaine network. The 
average (open circles) and the standard déviation (solid squares) are calculated 
from 100 square areas A, randomly chosen within the network. The average minus 
standard déviation (down triangles) shows no variation for areas between 1 and 
5000 km . The standard déviation follows a power law (solid line) for areas 
between 1 and 100 km 2 . 

Figure 14. Scale dependency of the fork density. Same symbols as for figure 13. 

Figure 15. Plot of the mainstream length versus drainage area for the Vilaine (top) 
and Blavet (bottom) networks. The best power-law fit is drawn with a solid line; 
the dashed line show the asymptotic behaviour (Z M ~ VÂ). 
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Table 1 



•Kb R] «s a z «a a M a H « T 

Vilaine 4.1-4.3 2.3-2.5 1.65 1.42 1.30 1.43 0.65 0.8 

1.6-1.7 1.40-1.44 1.25-1.35 1.41-1.45 



Blavet 4.0-4.2 2.1-2.4 1.75 1.45 1.35 1.7 0..65 0.85 

1.65-1.85 1.43-1.47 1.3-1.4 1.5-1.75 
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CHAP. 2: MODÈLES DE CROISSANCE DES RÉSEAUX 
HYDROGRAPHIQUES ET D'ÉVOLUTION DU RELIEF 



1. Introduction 

L'objet du chapitre 1 était de contraindre la dynamique des reliefs et des réseaux 
hydrographiques qui leur sont associés, par une analyse de cas de figures naturelles, afin 
de faire ressortir leurs caractéristiques spatiales propres. Cette analyse a l'inconvénient de 
faire implicitement référence à des échelles de temps et d'espace finies que sont la durée 
de vie et la taille du système étudié. Mais quand bien même l'étude d'objets réels aboutit à 
un schéma d'organisation, il est nécessaire d'une part, d'approfondir ce schéma en 
validant des hypothèses de mise en place des réseaux et de topographies à l'aide de 
modèles d'organisation, et d'autre part, d'utiliser ces derniers dans une démarche 
d'extrapolation du comportement des systèmes naturels sur des plages d'échelles de 
temps et d'espace plus étendues. Ceci est l'objet de ce chapitre. 

Si l'analyse de la structure géométrique des réseaux hydrographiques et des 
propositions de modèles d'organisation furent développées dès les années 40, la 
validation et l'utilisation de ces résultats furent contraintes par le développement des 
moyens informatiques. L'avènement de puissance de calcul et de mémoires numériques 
conséquentes a permis entre autres, l'utilisation et l'exploitation de topographies 
numériques, et par là même un renouveau dans l'étude quantitative et la modélisation des 
structures géomorphologiques. Nous pouvons distinguer deux grands groupes de 
modèles géomorphologiques proposés dans la littérature : 

(i) les modèles de réseaux qui s'attachent à reproduire des arborescences 
structurées similaires aux réseaux naturels, à l'aide soit de critères 
d'organisation des connections {modèles de configuration), soit de règles de 
croissance dans un espace bi-dimensionnel {modèles de croissance), 

(ii) les modèles couplés réseau-topographie simulant la formation de reliefs et de 
réseaux hydrographiques basés sur une modélisation des processus d'érosion. 

Le premier type de modèles nous permet de discuter les règles d'organisation 
observées dans le chapitre 1. Le deuxième type de modèles nous permet d'aborder la 
dynamique du système mécanique et le développement des instabilités les plus 
représentatives. Nous développerons, pour chacun de ces deux groupes, les modèles que 
nous estimons les plus représentatifs et nous proposerons en fin de chapitre un modèle 
d'érosion basé sur nos propres observations. 
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2. Les modèles de réseaux 

L'évolution des réseaux hydrographiques suit selon certains auteurs (Glock, 1931; 
Howard, 1990) différentes phases que l'on peut grossièrement réduire à : une première 
phase d'extension et d'occupation de l'espace disponible, limitée à une densité maximum 
de drainage en relation avec les conditions locales (pluviométrie, végétation, géologie), 
puis une seconde phase de réorganisation du réseau (modification de la topologie des 
branchements) pendant laquelle le système tend vers une configuration d'équilibre. 

Les modèles, présentés ci-dessous, diffèrent par l'importance donnée à la phase de 
croissance sur la configuration d'équilibre. Les premiers modèles (modèles 
d'optimisation) considèrent une croissance aléatoire et une topologie des branchements 
fortement modifiée pendant la phase de réorganisation. Les autres modèles considèrent 
que l'organisation est atteinte à la fin de la phase de croissance. 

Les modèles suivant tentent, à partir de règles déterministes ou stochastiques, de 
retrouver les structures des réseaux hydrographiques, et proposent ainsi des principes 
d'organisation possibles. 

2.1. Les modèles d'optimisation 

Sous cet intitulé nous regroupons des modèles de construction de réseaux 
hydrographiques basés sur un principe d'optimisation. Ces modèles définissent un critère 
à minimiser et ajustent la structure du réseau tant que la configuration optimale n'est pas 
atteinte. Deux modèles à notre connaissance dominent dans ce groupe: le modèle 
d'Howard (1990) basé sur une minimaUsation de la puissance par unité de longueur de 
cours d'eau, et le modèle de Rodriguez Iturbe et al. (1992ab) basé sur une minimalisation 
locale et globale de l'énergie dissipée. 

2.1.1. Géomorphologie et thermodynamique 

Un paysage peut être considéré comme un système ouvert soumis à l'action 
antagoniste de processus endogènes dus aux mouvements tectoniques, et des processus 
exogènes comme la pluie ou le vent. Les processus endogènes apportent de la matière aux 
systèmes (flux de matière > 0), au contraire des processus a<ogènes dont l'action tend à 
retirer de la matière au système (flux < 0). L'action simultanée de ces deux types de 
processus conduit-elle vers un état d'équilibre du système? La détermination des 
conditions d'équilibre, et l'influence de perturbations sur la stabilité des états d'équilibre, 
constitue l'enjeu principal de la compréhension de la dynamique des reliefs. Cette 
approche est très similaire à celle de la thermodynamique classique ou statistique. Il est 
donc assez naturel de tenter d'effectuer un rapprochement, voire une assimilation, entre 
les grandeurs thermodynamiques et géomorphologiques, et d'établir les mêmes lois de 
comportement. 
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L'analogie avec un système thermodynamique fut introduite par Leopold et Langbein 
(1962), en identifiant les variables géomorphologiques à des variables 
thermodynamiques. Ils proposèrent ainsi la correspondance suivante: 

T<->h 
dQ<->dM 

où h, l'altitude par rapport à un niveau de base donné, joue le rôle d'une température 
thermodynamique, et dM, la variation de quantité de masse, est analogue à dQ, la 
variation de quantité de chaleur. Nous pouvons dès lors définir l'entropie S d'un système 
thermodynamique ou d'un relief comme: 

dS=dQ/TodM/h. 

L'intérêt de l'analogie thermodynamique est de définir l'état d'équilibre du système à 
partir des principes de la thermodynamique dont le deuxième est que le taux de 
production d'entropie est minimal à l'équilibre. 

L'application de ces principes thermodynamiques généraux permet ainsi de retrouver 
le profil en long caractéristique des rivières (h ~e- cx , avec x la distance par rapport à la 
ligne de crête et c une constante), et donne ainsi une certaine crédibilité à cette approche 
(Scheidegger et Langbein, 1966; Scheidegger, 1967). Une autre approche, très similaire 
à la précédente, consistera à minimiser l'énergie dissipée, afin de trouver l'état du 
système le plus stable. L'effort de modélisation portera dans ce cas non pas sur la 
détermination de l'entropie, mais plutôt sur l'expression de l'énergie dissipée par les 
processus d'érosion. 

2.1.2. Modèle d'optimisation des angles de confluence (1990) 

D'après Howard quatre processus importants modifient l'aire drainée des bassins 
versants: 

(1) le développement de nouveaux biefs par une croissance des ordres 1, 

(2) la variation de l'aire drainée par une migration des lignes de crête, 

(3) les processus de capture, 

(4) les migrations de vallées et les réajustements d'angles de confluence. 

Si le processus (1) intervient très tôt dans l'élaboration des réseaux hydrographiques, 
les processus (2) (3) et (4) opèrent de façon beaucoup plus lente. Les corrélations 
observées sur des réseaux hydrographiques entre les angles de jonction des biefs avec le 
débit Q et la pente, laissent suggérer l'existence d'une configuration optimale des 
écoulements au niveau des confluences. Une érosion différentielle des berges 
engendrerait une migration des vallées et par la-même une évolution des angles initiaux de 
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confluences vers une configuration plus optimale. Le modèle d'Howard s'appuie sur 
l'idée que le taux de travail des écoulements Q, au niveau de chaque confluence, doit être 
minimal. 

3 

i=l 

avec i le nombre de biefs, Cj la puissance par unité de longueur de bief, L; la longueur 
du bief i. En posant: 

Q = pgQiSi 
et 

Si = KQf = K{aLif 

avec Si est la pente moyenne du bief i, p la densité du fluide, K est une constante, OC le 
taux d'apport en eau des berges (supposé constant sur l'ensemble du réseau), Q ne 
dépend plus que de deux variables à savoir la longueur des biefs et le paramètre z. La 
structure optimale est obtenue par un calcul itératif des valeurs de 13 . A chaque pas du 
calcul le réseau est réactualisé par un déplacement des confluences d'une quantité 
proportionnelle à la longueur moyenne des biefs et dans le sens opposé à VX2. 
L'application de ce modèle à un réseau théorique obtenu par une procédure de croissance 
aléatoire, permet d'obtenir un réseau d'aspect plus réaliste (Fig II 1). Notons cependant 
que l'aspect du réseau final dépend fortement de z: lorsque z — » 0 le réseau présente des 
directions de drainage plus alignées avec des angles de confluence plus faibles. Les 
réseaux résultats présentent cependant les caractéristiques topologiques et géométriques 
des réseaux naturels. 

L'une des faiblesses de ce modèle, comme le reconnaît son auteur, est d'attribuer une 
valeur constante au paramètre z sur l'ensemble du réseau et du calcul d'optimisation. 
L'expression de la pente en fonction du débit devrait en effet dépendre de l'échelle du 
bassin versant étudié et évoluer au cours de la migration des vallées. Le modèle peut donc 
être optimisé en changeant la valeur de z suivant la valeur de l'aire drainée amont. Mais 
quand bien même l'optimisation de ce modèle soit réalisable, il reste à trouver d'autres 
critères de validation permettant de faire ressortir l'efficacité de l'optimisation. Le modèle 
initial vérifie en effet déjà les critères de validation couramment utilisés. 
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L j is the total length of the paths, 
and Lis the average length of the path from a point to the outlet. In 
the explosion case, L'f ' refers to the case when there is a minimum 
displacement among the points so that„tbsre js a différent path 
between each point and the outlet ,' 



Figure II 2 : 3 types de schémas de connections possibles. Si l'énergie dissipée par le système est 
proportionnelle à la longueur du chemin de connexion, le type arborescent est la structure qui dissipe 
à la fois localement et globalement le moins d'énergie, (d'après Stevens, 1974) 






Figure II 3 : Optimal Channel Network (OCN) obtenus à partir d'une configuration aléatoire a). 
L'aire drainée critique nécessaire à la formation d'un bief est dans cet exemple égale à 4 pixels. 

Critère d'optimisation: J Pj = % 2 j^L,-. 

i i 

Le réseau b) correspond au cas h=0.5, le cas c) à h=l. (d'après Rigon et al, 1993). 



2.1.3. Modèle des "réseaux optimisés" ou "Optimal Channel Networks" (OCN) 
(1992-) 

Le réseau hydrographique est ici considéré comme un système dissipatif qui d'après 
les lois thermodynamiques doit tendre vers un état où l'énergie dissipée est rrrinimale. Si 
l'on suppose que l'énergie dissipée P dans un bief dépend linéairement de sa longueur, on 
démontre que le système arborescent est celui qui dissipe localement et globalement le 
moins d'énergie (Stevens, 1974; Fig. II 2). 

Inspirés par ce principe, Rodriguez Iturbe et al. (1992a) ont construit un modèle basé 
sur trois postulats: 

(1) le taux d'énergie dissipée par chaque bief doit être minimal; 

(2) le taux d'énergie dissipée par unité d'aire est constant sur toute la surface du 
bassin versant; 

(3) le taux d'énergie dissipée totale du système doit être minimale. 

L'application des deux premiers postulats conduit à un système de drainage 
arborescent simple, où la vitesse d'écoulement est constante sur tout le réseau 
hydrographique. L'énergie dissipée peut alors s'exprimer en fonction du débit moyen 
annuel Qi et de la longueur L/ de chaque bief: Pi = kQ°i' 5 Li. La condition de minimisation 
imposée par le troisième postulat peut alors s'écrire: 

E =XP/ = X^°X- 

Le débit Qi sera par la suite assimilé à l'aire drainée Ai dans les calculs. La simulation 
proprement dite consiste à chaque pas d'itération, à choisir un lien au hasard, à changer 
son orientation, et à retenir pour l'itération suivante la configuration pour laquelle E prend 
la plus faible valeur. L'expression de l'énergie E est similaire à celle définie dans le 
modèle d'Howard. Les deux modèles se distinguent essentiellement par l'application 
d'une règle de minimalisation de l'énergie dissipée sur des échelles d'espace différentes. 

Les réseaux théoriques obtenus vérifient l'ensemble des critères topologiques et 
géométriques des réseaux naturels (Fig. II 3 ; Rinaldo et al., 1992; Rigon et al., 1993). 
L'expression de l'énergie dissipée E reste néanmoins l'un des points sensibles du 
modèle. L'utilisation d'une expression plus générale: E = ^ kAfLi , montre que les 

propriétés topologiques des réseaux résultats sont préservées si 0<n<l (Rinaldo et al., 
1992). Il serait par conséquent, comme dans le modèle précédent, intéressant d'étudier 
l'effet d'une modification de l'hypothèse simplificatrice à l'origine de l'expression 
E = ^ kQl' 5 Li , et selon laquelle la vitesse des écoulements est supposée constante sur 

i 

tout le réseau hydrographique. 
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Les OCN ont la particularité: 

(i) de ne pas dépendre des conditions initiales, 

(ii) d'aboutir à des résultats semblables lorsque le système est soumis à des 
perturbations; 

(iii) de présenter une distribution des aires drainées suivant une loi de puissance 
(n(A) = A 0 ' 410,1 ) caractéristique des systèmes d'agrégation décrits par 
Takayasu et al. (1988) et observée sur des cas naturels (chapitre 1). 

Ces trois propriétés des OCN démontrent selon les auteurs de ce modèle, qu'un principe 
de minimalisation de l'énergie dissipée appliqué à la fois à une échelle locale et globale 
peut engendrer des systèmes dissipatif tendant naturellement vers un état critique (Rinaldo 
et al., 1992). L'état critique d'organisation correspond à une répartition des aires drainées 
suivant une loi de puissance sans longueur caractéristique. Autrement dit, le système 
s'organise de lui-même, quel que soit l'état initial et quelles que soient les perturbations 
auxquelles il est soumis, vers une configuration géométrique telle que la distribution des 
aires drainées suit une loi de puissance avec une puissance égale à 0,4 ±0,1. Cet état 
critique n'est pas synonyme d'état d'équilibre, mais simplement d'une configuration type 
vers laquelle le système tend et conserve au cours de son évolution. D'après Rinaldo et 
al. (1992) cette dynamique d'évolution serait à l'origine des propriétés fractales des 
réseaux hydrographiques, ces dernières représenteraient en fait les principaux indices en 
faveur de ce modèle. 

Poursuivant ce raisonnement Rinaldo et al. (1993) et Rigon et al. (1994) proposent un 
modèle d'érosion inspiré de la dynamique particulière des tas de sable. Un tas de sable est 
en effet l'exemple type de système auto-organisé tendant vers un état critique (Bak et al. 
1987,1988; Christensen et al., 1992). A partir d'une configuration d'altitude et d'un 
réseau de drainage initiaux, on calcule en chaque point de la topographie la contrainte 
locale cisaillante t,- engendrée par un écoulement de surface de débit local Qf. 

Ti = pgQ?> 5 iAz =pgA?' 5 iAz 

avec p est la densité de l'eau, g l'accélération liée à la gravité, Az la différence d'altitude 
suivant la ligne de plus grande pente. Le principe est de considérer une contrainte 
cisaillante seuil Te au- delà de laquelle le système est localement instable (déclenchement 
d'une avalanche). A chaque pas d'itération, le pixel présentant une différence (t,- - Te ) 
maximale voit son altitude diminuée de manière à imposer t,- = Te . La matière érodée est 
retirée du système sans sédimentation. Le réseau de drainage est alors recalculé pour 
établir la nouvelle grille de contraintes t,- . On répète cette procédure tant qu'il existe un 
point tel que t > Te . Le système évolue ainsi vers un état "critique". Puis on introduit 
une perturbation en modifiant l'altitude d'un point choisi au hasard. Cette perturbation 
induit une nouvelle configuration du réseau de drainage et une nouvelle grille des 
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contraintes ce qui engendre une nouvelle phase d'avalanches en cascade. On perturbe 
ainsi le système tant que la perturbation produit une modification du réseau de drainage. 
Le réseau de drainage résultat vérifie les propriétés topologiques et géométriques des 
réseaux hydrographiques, et valide cette approche selon les auteurs de ce modèle. 

Les modèles d'optimisation ont le mérite de montrer qu'une approche 
thermodynamique conduit à des réseaux réalistes. Cependant, l'hypothèse selon laquelle 
les réseaux hydrographiques évoluent d'un état initial vers une configuration optimale, 
génèrent des systèmes insensibles aux conditions initiales et très stables vis à vis 
d'éventuelles perturbations. Ces propriétés, contradictoires avec les résultats obtenus par 
des simulations analogiques (Schumm et al., 1987), représentent un des points faibles de 
ces modèles. La condition de minimalisation de l'énergie dissipée apparaît donc comme 
une contrainte trop forte dans l'évolution d'une topographie. 

2.1.4. Modèle topologique aléatoire (1966) 

Les modèles aléatoires, à travers une vision stochastique des processus, constituent un 
pôle important dans l'ensemble des modèles de construction de réseaux hydrographiques. 
Le terme aléatoire peut être ici perçu sous deux angles différents: 

(i) les processus étudiés comportent une composante aléatoire intrinsèque qui les 
rendent imprévisibles à plus ou moins long terme; 

(ii) la complexité ou le nombre des processus impliqués ne permet pas à 
l'observateur, avec les moyens matériels dont il dispose, de saisir tout le 
déterminisme du système. 

Partant de la simple constatation que le nombre de sources n d'un réseau 
hydrographique fixe le nombre de confluences (n-1) et le nombre de liens entre deux 
confluences (2n-l), Shreve (1966) émet l'hypothèse que des réseaux de même magnitude 
(même nombre de sources) sont d'égale complexité avec une probabilité identique 
d'exister. Dans le modèle de Shreve, tout réseau naturel peut être considéré comme l'une 
des réalisations possibles d'une population de réseaux de topologie aléatoire (le terme 
topologie fait ici référence au schéma de connexions utilisé pour relier les confluences 
entre elles, la distance entre deux confluences est supposée constante). Une telle 
population regroupe tous les réseaux de même magnitude et définit un ensemble de 
réseaux topologiquement distincts (Topologically Distinct Channel Networks: IDCN ; 
Fig. II 4). Le nombre de TDCN de magnitude n au sein d'une même population est 
obtenu par la formule: 

2n-\ ni {n-ly. 
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Par une analyse combinatoire des populations de TDCN, Shreve (1966) démontre que 
la configuration la plus probable vérifie les lois de Horton (Rb=4, Rl=2). 

Ce modèle sera par la suite généralisé (Shreve, 1967, 1969, 1975; Smart, 1968; Smart 
et Wemer, 1976) par deux postulats: (i) en l'absence de contraintes d'environnement les 
réseaux hydrographiques ont une topologie aléatoire; (ii) les biefs intérieurs et extérieurs 
et les aires drainées qui leur sont associées ont des distributions statistiques 
indépendantes, qui ne dépendent pas de leur localisation au sein du réseau. 

Le modèle de Shreve est l'un des modèles qui concordent le mieux avec les 
observations. Les résultats statistiques obtenus sur l'analyse de la Vilaine et du Blavet 
confirment l'hypothèse d'une configuration topologique et géométrique aléatoire 
(distribution exponentielle des longueurs de biefs intérieurs et extérieurs). 

Ce modèle n'est pas à proprement parlé un modèle d'optimisation car il présuppose 
l'absence de déterminisme dans la construction du réseau. Cette hypothèse, quoique très 
insatisfaisante sur le plan de l'analyse des processus, met néanmoins en évidence de façon 
théorique l'éventualité d'une croissance aléatoire des réseaux hydrographiques. 

2.2. Les modèles de croissance 

La deuxième grande catégorie de modèle part implicitement de l'hypothèse que la 
phase de croissance des réseaux fixe la topologie des structures hydrographiques, n 
existe un certain nombre de modèles qui, à partir de cette hypothèse, proposent des règles 
d'extension des arborescences hydrographiques indépendamment des processus 
d'érosion, avec deux approches complémentaires : 

(i) l'extrapolation d'observations effectuées sur des petits bassins versants ( < 1 
km 2 ), ou des modèles réduits analytiques, aux échelles supérieures, 

(ii) l'analogie avec des processus de croissance de structures similaires observés 
pour des systèmes de nature différente. 

L'étude de la croissance de réseaux de taille restreinte conduit à des modèles 
déterministes très qualitatifs présentant un développement en étapes distinctes. Difficiles à 
transcrire sous forme numérique, ces modèles ont néanmoins le mérite de décrire en 
quelques traits principaux le développement de réseaux expérimentaux, et trouvent à ce 
titre leur place dans ce chapitre. L'assimilation des réseaux hydrographiques à des 
systèmes d'agrégation, oriente la recherche de modèle de croissance vers des modèles 
d'agrégation de molécules. Nous présenterons dans ce texte deux modèles d'agrégation 
qui permettent de reproduire des arborescences satisfaisant aux critères de validation 
courants, à savoir l'invasion percolation et l'agrégation limitée par la diffusion. 
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2.2.1. Les modèles déterministes 



2.2.1.1. Modèle de Horton (1945) 

Le modèle conceptuel de croissance de Horton s'attache à décrire l'extension des 
réseaux hydrographique dans le contexte particulier d'une surface uniforme soumise à un 
ruissellement généralisé. Horton suggère qu'au-delà d'une certaine distance critique xc, 
des rigoles parallèles au sens de la pente apparaissent et commencent à inciser la surface. 
Les rigoles dont le taux d'incision est le plus fort vont augmenter latéralement l'aire 
qu'elles drainent, capturer les rigoles voisines et initier ainsi un réseau dendritique (Fig. 
II 5a). 

Ce modèle très simple ne permet pas de décrire la formation des réseaux 
hydrographiques sur des surfaces présentant des géométries et des propriétés variées 
(rugosité, végétation, perméabilité). Il se limite au cas de figure d'une surface plane peu 
pentée et présentant une perméabilité faible. De plus Horton ne donne pas d'explication 
sur l'origine physique de l'apparition du réseau parallèle de rigoles. L'utilisation d'un 
seuil de comportement (longueur critique xc), constitue le point le plus intéressant du 
modèle avec l'introduction d'une vision non linéaire des processus mis enjeu. 

2.2.1.2. Modèle d' érosion progressive (1956) 

Schumm (1956), Howard (1971), Smart et Moruzzi (1971) proposent un modèle où le 
réseau hydrographique progresse d'une zone drainée par un réseau dense vers une zone 
vierge de tout réseau (Fig. II 5b). La progression s'effectue par un allongement et une 
subdivision des ordres 1 qui, pris dans leur ensemble, peuvent être ainsi assimilés à une 
onde de dissection de la surface disponible. 

Ce modèle, comme le précédent, reste très qualitatif. Sa particularité est de supposer 
une progression purement aval-amont du réseau hydrographique, sans partage des 
surfaces à l'amont des biefs d'ordre 1. 

2.2.1.3. Modèle de Glock (1931 ) 

Glock propose un modèle basé sur une croissance séquentielle du réseau 
hydrographique. Il distingue 4 phases d'évolution (Fig. II 5c): 

(1) Initiation: le réseau se présente sous la forme d'une distribution éparse et 
discontinue de biefs; 

(2) Extension: élongation et raccord des biefs initiaux Ç'squelette" du réseau), 

(3) Maximum extension: ramification sensible de l'arborescence de base afin d'occuper 
tout l'espace disponible; 

(4) Intégration: élimination des biefs dont les bassins versants, par aplanissement du 
relief, ne peuvent plus entretenir l'existence; 
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Figure II Ç : Modèles conceptuels déterministes d'évolution de réseau hydrographique. Le modèle A 
correspond au modèle de Horton: un réseau parallèle de rigoles est progressivement remplacé par 
une structure arborescente; le modèle B correspond à une progression des biefs d'ordre 1 vers la 
partie amont du bassin non encore drainée; Le modèle C considère une installation rapide d'un réseau 
principal complété dans une deuxième phase par la formation de branches secondaires, (d'après 
Schumm et al., 1987) 
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Figure II 4 : Les 14 Topologically Dinstinct Channel Networks (TDCN) de magnitude 5 (nombre de 
sources = 5) (d'après Shreve, 1966) 
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directions 



Marcheur aléatoire 



jfjgjj maille de réseau 
> parcours du marcheur 



Pi probabilité d'aller dans la direction i 



Figure H6 ; Technique du marcheur aléatoire (exemple de la méthode D4: 4 directions possibles de 
déplacement). Un marcheur lancé sur une topographie. Le chemin suivi par le marcheur est 
conditionné par les probabilités Pi de déplacement. Les valeurs de Pi dépendent des conditions 
locales (pente, perméabilité...). Si Pi est constante quelque soit i, le parcours est aléatoire. 



Ces quatre étapes d'évolution correspondent à une description trop qualitative des 
processus pour pouvoir être appliquée de façon reproductible à des cas concrets. Les 
observations sur des réseaux réels ne permettent pas en effet de distinguer clairement ces 
différentes étapes. 

Les trois modèles déterministes précédents constituent des modes différents 
d'extension des réseaux hydrographiques. Si l'on se réfère aux expériences analogiques 
effectuées au Colorado State University (pour plus de détails, le lecteur pourra se 
référencer à Schumm et al., 1987), le modèle de Glock apparaît comme le mode 
d'évolution le plus courant. Un des faits majeurs mis en évidence lors de ces expériences 
est l'importance des conditions initiales. La structure finale du réseau hydrographique et le 
mode de croissance semblent en effet très dépendants de paramètres comme la valeur 
initiale de la pente ou du niveau base. De manière générale une variation minime des 
conditions initiales engendre des modifications importantes dans la structure des réseaux 
développés. 

2.2.2. La marche aléatoire 

2.2.2.1. La technique du marcheur 

Une autre grande famille de modèles utilise la technique du marcheur pour générer des 
réseaux théoriques. Le parcours d'un marcheur dans un espace comme une topographie 
ou un système maillé, peut être conditionné par des règles plus ou moins contraignantes. 
Ces dernières définissent, à chaque "pas" d'itération et pour chaque direction possible, 
une probabilité de déplacement (Fig. II 6). L'expression de la probabilité de déplacement 
en fonction de paramètres physiques introduit dans le modèle un déterminisme 
correspondant à l'application de règles physiques ou stochastiques vis à vis des propriétés 
locales du système. Les caractéristiques physiques du système (valeur d'un paramètre 
pour chaque unité du système) se présentent sous la forme de "topographie" propre à 
chaque paramètre, et sur lesquelles il est possible de calculer des potentiels et par la-même 
des gradients qui vont conditionner la marche. Parallèlement aux règles de marche, des 
règles d'action du marcheur sur la ou les "topographies" peuvent être définies 
(modification de la surface sur laquelle se déplace le marcheur). Le "parachutage" des 
marcheurs peut s'effectuer de façon ordonnée dans le temps ou dans l'espace, ou à 
l'inverse être complètement aléatoire. 

La technique du marcheur est un moyen très souple d'étude de comportement des 
systèmes spatialisés. Au-delà de sa forme assez intuitive (on identifie aisément le 
marcheur à une goutte ou à une certaine quantité d'eau), le grand avantage de cette 
méthode est de pouvoir introduire simplement une composante aléatoire plus ou moins 
forte dans les modèles. Elle demande cependant une certaine réflexion dans 
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l'interprétation des équations déterministes qui gouvernent l'évolution du système, en 
termes de règles de conduite du marcheur. 

2 .2 .2 .2 . Modèle du marcheur ivre (1962 ) 

Les modèles de marcheurs utilisés pour construire des réseaux de drainage dans un 
plan, ont pour objectif de définir les connections reliant les nœuds d'un système maillé 
symbolisant les confluences ou les points sources d'un réseau hydrographique. Précisons 
que cette procédure qui suit ne permet pas d'étudier les processus de croissance, et simule 
simplement des configurations possibles d'arborescences. On impose que des règles de 
marche, la goutte d'eau ne modifiant pas sur son passage les propriétés de la surface de 
déplacement. Lorsqu'un marcheur est lancé sur la grille, il poursuit sa course tant qu'il 
n'a pas atteint les limites du système ou un site déjà visité. Le chemin suivi par chaque 
marcheur représente un des biefs du réseau final. Afin d'éviter la formation de réseaux 
bouclés, on empêche le marcheur de recouper sa course (Self-Avoiding Walk: SAW). La 
marche est ici aléatoire: à chaque itération, la probabilité de déplacement du marcheur est 
égale quelles que soient les directions possibles. La figure II 7 représente l'évolution du 
réseau après " le lancé" de 3 marcheurs. 

Leopold et Langbein (1962) proposèrent les premières simulations de SAW pour 
reproduire la structure des réseaux hydrographiques (Fig. II 8). Les réseaux théoriques 
obtenus respectent les lois de Horton et vérifient approximativement la loi de Hack. 
Meakin et al. (1991) étudièrent la sensibilité des SAW vis à vis des règles de parcours. Ils 
proposent notamment une autre variété de marches aléatoires, où les sites cibles du 
marcheur ne sont plus limités aux plus proches voisins, mais choisis suivant une direction 
aléatoire (Indefinetely Self-Avoiding Walk: IGSAW). Quel que soit le type de trajectoire, 
les lois de Horton et la loi de Hack restent satisfaites avec une marge d'erreur raisonnable. 
Les résultats obtenus dans le cas des SAW mettent en évidence des effets de taille finie, et 
montrent un comportement asymptotique des lois géométriques avec la taille du système. 
Dans tous les cas, les réseaux générés par une marche aléatoire vérifient la loi de 
puissance propre au modèle d'agrégation de Takayasu et al. (1988). Notons pour 
mémoire que Scheidegger (1967) et Seginer (1969) présentent également des SAW très 
semblables de part leur construction et leurs résultats aux deux modèles précédents. 

Les modèles stochastiques reproduisent les lois géométriques et statistiques observées 
pour les réseaux naturels (lois de Horton, loi de Hack, loi de puissance des aires 
drainées). Es démontrent ainsi qu'il n'est pas nécessaire d'imposer un critère de 
minimisation pour satisfaire les lois d'organisation observées. On peut donc s'interroger 
sur la représentativité des critères de validation utilisés. Les travaux de Kirchner (1993, 
1994a, 1994b) mettent en évidence l'aspect très général des lois de Horton, en démontrant 
qu'une arborescence binaire quelconque (deux biefs par confluence) vérifie ces lois. De 
même, le modèle d'agrégation de Takayasu, sur lequel s'appuient les conclusions de 
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a) 




Figure II J : Construction d'un réseau théorique à l'aide de la technique du marcheur aléatoire. Les 
figures a), b) et c) représentent 3 stades de construction pour 3 lancés de marcheur successifs. Le 
parcours d'un marcheur aléatoire constitue la formation d'une branche du réseau théorique. 




Figure II %: Réseau théorique obtenu à l'aide de marcheurs aléatoires 
(d'après Leopold et Langbein, 1962). 



Rinaldo et al., n'utilise aucune règle déterministe d'agrégation et de ce fait considère ce 
phénomène de façon globale. C'est pourquoi, il est prématuré d'une part, d'attribuer à ces 
critères un caractère propre aux réseaux hydrographiques, et d'autre part, d'interpréter ces 
tendances comme la signature d'un phénomène déterministe. L'impossibilité de pouvoir 
déterminer l'efficacité relative des modèles à l'aide de critères quantitatifs ne permet pas de 
conclure sur la nature plus ou moins aléatoire des processus. Les modèles d'organisation 
aléatoires s'avèrent cependant les modèles les plus simples satisfaisant les critères de 
validation actuels, mais également ceux dont les principes s'accordent le mieux avec les 
lois statistiques observées (Chapitre I). L'hypothèse d'une réorganisation du réseau 
hydrographique suivant des règles déterministes après une première phase d'extension, ne 
présente pas d'argument supplémentaire par rapport à l'hypothèse d'une construction 
aléatoire. La structure aléatoire des longueurs de cours d'eau sur des échelles supérieures 
à 1 km, laisse au contraire envisager la construction d'une arborescence stable pendant la 
phase de croissance. 

2.2.3. La percolation d'invasion appliquée aux réseaux hydrographiques 
2.2.3.1. La percolation d invasion 

L'invasion percolation est à l'origine un modèle mis au point pour simuler la 
progression de fluides dans un milieu poreux (Meakin, 1991). Cette méthode consiste 
dans un premier temps, à définir sur un système maillé un champ aléatoire de cohésion du 
système (porosité, résistance, composition chimique) par rapport à une contrainte 
extérieure (pression d'un fluide, force, corrosion). La valeur locale du champ, i.e. sur 
chaque maille, permet de définir une probabilité locale de réaction du système vis à vis du 
processus étudié. L'invasion du système s'effectue à partir d'un site ou d'un ensemble de 
sites (front de percolation) définis préalablement, qui déterminent la frontière entre les 
éléments du système qui sont inchangés et ceux modifiés sous l'action de la contrainte 
(pores vides/pores remplis d'eau, sites intacts/ sites cassés, site propres/sites rouillés). La 
règle de progression du front est la suivante: à chaque pas d'itération, le site en contact 
direct avec le front de percolation, possédant la probabilité la plus faible, est envahi et 
devient partie intégrante du front. La figure II 9 illustre trois itérations successives de 
percolation d'invasion. 

Une des propriétés des systèmes de percolation d'invasion est de générer des 
structures de dimension fractale D =1.82. Cette valeur, très proche de celle obtenue par La 
Barbera et Rosso (1989,1990) en utilisant les lois de Horion, encouragea l'application de 
cette méthode pour les simulations de réseaux hydrographiques. 
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Figure IL9 : Règle d'invasion percolation. Une structure initiale croît par occupation du site mitoyen 
possédant la plus faible résistance. L'invasion percolation présentée dans l'exemple s'effectue suivant 
8 directions. 
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Figure H. 10 : 

Examples of simulated drainage networks grown by the self-avoidir.,; 
invasion percolation method. Numerical simulations were set up using a 
square site-bond lattice in which random bond strengths were normaliy 
distributed. Lattice size is 512 x256. Seeding was allowed from svery other 
point along the bottom edge. In this case invasion of the lattice was allowed 
to continue until no further free sites were available. Sirahler 28 stream 
ordering was appliedt streams of order four and above are shown. The 
substrate becomes steadily more susceptible to northward growth from 3 
(isotropie lattice) to c (strongly anisotropic lattice). The fractal dimension 
o* min of the principal streams is invariant to this anisotropy and is on average 
1.130 throughout. 



(d'après Stark, 1991Ï 



22.3.2. Le modèle de Stark (1991) 

Stark (1991) propose d'appliquer un modèle d'invasion percolation pour simuler la 
croissance des réseaux hydrographiques. H apporte cependant une modification au modèle 
de manière à empêcher la formation de boucles: un site ne peut pas être envahi si l'un des 
ces voisins, excepté le site source, est déjà occupé (self-avoiding invasion percolation). 

Les réseaux théoriques (Fig. II 10) résultants vérifient les lois de Horton et prédisent 
avec une remarquable précision la loi de Hack (L~A a > ^simulation" 

0.565, 

ûexpérimental=0.568). Notons cependant qu'une telle méthode ne peut pas créer des 
réseaux de dimension égale à 2, et ne peut donc pas reproduire les résultats obtenus pour 
la Vilaine et le Blavet 

La percolation d'invasion peut être assimilée à une croissance des réseaux 
hydrographiques gouvernée par un contrôle aval. Si cette image correspond assez bien au 
remplissage d'un lac, elle reste néanmoins trop simpliste pour simuler l'action érosive des 
écoulements amont-aval présents dans les systèmes naturels. 

2.2.4. L'agrégation limitée par la diffusion appliquée aux réseaux 
hydrographiques 

22.4.1. L'agrégation limitée par la diffusion 

Utilisée pour simuler les processus d'agrégation irréversibles entre particules, 
l'agrégation limitée par la diffusion (Diffusion-Lirnited Agrégation: DLA) utilise toutes les 
possibilités offertes par la technique du marcheur (Meakin, 1991; Wicseck, 1992). 

Considérons une particule immobile sur la grille représentant le champ de déplacement 
des marcheurs. Cette particule sera une cible. Au-delà d'une distance r de la cible 
représentant l'infini, nous lançons un marcheur que nous laissons errer au hasard sur la 
grille mais avec une restriction essentielle: si le marcheur passe sur l'une des cases 
voisines de la cible, elle colle à la cible. On construit ainsi une nouvelle cible à deux 
particules. Puis nous lançons un autre marcheur... La figure II 1 1 présente l'évolution du 
système après le lancé de 3 marcheurs. 

A l'instar des figures de percolation, les agrégats produits par un DLA possèdent une 
dimension bien définie Z)=l,65. Cette dimension correspondant au cas précédent, peut 
varier suivant les règles de parcours et de "collage". La variété de règles possibles offre 
un potentiel infini de cas de figures, et autorise la simulation numérique d'un grand 
nombre de processus de transfert. 
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maille de réseau 



parcours du marcheur 



Figure II. 1 1 : Construction d'un réseaux théorique par agrégation limitée par la diffusion. Les 
schémas a), b) et c) représentent 3 stades successifs de croissance. Le réseau se développe par 
agrégation de marcheurs aléatoires sur une structure initiale qui peut être un point ou une 
arborescence. 




Figure 11.12 : Exemple de réseaux théoriques construits par agrégation limitée par la diffusion. Dans 
cet exemple on impose une contrainte locale d'agrégation empêchant les marcheurs de "coller" sur 
des sites adjacents à deux cellules du réseau, (d'après Masek et Turcotte, 1993) 



2.2.42. Application aux réseaux hydrographiques (1993) 

Afin de mieux simuler la pluie et l'influence mutuelle exercée par la croissance des 
biefs (compétition entre les sources pour la partage de l'espace disponible), Masek et 
Turcotte (1993) présentent une version modifiée de DLA (Fig. II 12). Dans ce modèle: (i) 
les marcheurs sont lancés de manière aléatoire sur toute la grille et non plus à "l'infini", 
(ii) le collage du marcheur n'a lieu que dans les mailles ne possédant qu'un seul voisin 
appartenant à l'agrégat, (iii) le marcheur termine sa course s'il atteint les limites du 
système ou passe sur un site jouxtant l'agrégat. Les réseaux construits à l'aide de cette 
procédure vérifient les lois de Horton (R&=3,98; /?/=2,09) et possèdent une dimension 
fractale D=l,87. 

Bien que ce dernier résultat ne soit pas en accord avec les données sur la Vilaine et le 
Blavet, la technique du DLA garde un attrait certain de part sa souplesse d'utilisation et les 
possibilités d'amélioration qu'elle peut offrir. C'est pourquoi nous l'utiliserons par la 
suite comme schéma de base pour l'élaboration d'un modèle couplé réseau-topographie. 

2.3. Conclusion des modèles de réseaux 

Dans le chapitre 1, nous avons montré que les réseaux hydrographiques pouvaient se 
décrire par (i) une densité d'aire drainée constante, (ii) une distribution aléatoire des 
branchements à l'intérieur du réseau, (iii) des lois de puissance observées sur les 
distributions des paramètres qui intègrent toute l'information "amont" (aires drainées, 
longueur totale en amont...). L'ensemble de ces données peut s'interpréter simplement 
par un modèle de branchement aléatoire (Shreve, 1966) sur un arbre convergent 
(Takayasu et al., 1988). 

En fait, on peut montrer que la condition de branchement vers l'aval est suffisamment 
forte pour déduire l'ensemble des lois topologiques actuellement mesurées (lois de 
Horton, loi de Hack, distribution de longueurs...). Il n'est donc pas étonnant que les 
modèles, présentés précédemment, satisfassent aux lois topologiques des réseaux 
naturels. La faiblesse des études précédentes ne vient pas des modèles mais de l'absence 
d'identification de paramètres mesurables pertinents, qui ne soient pas la conséquence de 
la condition de branchement "vers l'aval". 

Les modèles d'agrégation, type DLA, prédisent des lois d'échelle spatiales qui ne sont 
pas observées sur les réseaux du Blavet et de la Vilaine décrits au chapitre 1. Cette 
contradiction peut avoir deux causes : soit l'analogie n'est pas pertinente et il faut 
invoquer d'autres mécanismes de croissance, soit il existe une limite à la phase de 
croissance "fractale" lorsque le réseau remplit tout l'espace disponible. En tout état de 
cause, il n'est pas possible de valider ce type d'analogie sans observation de réseaux non 
permanents et/ou sans modélisation directe incluant les processus d'érosion- 
sédimentation. 
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Figure H.13 : Évolution théorique de la pente d'un versant d'après la relation- 

-, dt 
a) F = constante; b) F = y; c) F = g. (d'après Scheidegger, 1961, 1991) 



3. Les modèles d'érosion 

La complexité des processus géomorphologiques vient de la superposition de 
mécanismes lents (altération chimique, transfert en milieux poreux, tectonique...) et de 
mécanismes rapides (transport dans les rivières, avalanches...). 

Les modifications du relief sont le résultat d'une chaîne de processus qui ne peuvent 
être correctement décrits que par des couplages fortement non-linéaires. 

Dans ce chapitre, nous décrirons les principales règles d'érosion, les mécanismes 
couplés qu'il faudrait prendre en compte pour une description réaliste de la dynamique de 
ces systèmes et les instabilités mécaniques qui en résultent 

3.1. Les équations générales d'évolution 

L'expression générale de l'évolution de l'altitude z d'un point (x,y) soumis à une 
action endogène tectonique et à une action exogène érosive, se déduit de l'équation de 
continuité de la matière en ce point: 



avec: V ,q s la divergence du flux de matière, et T(x,y,f) la composante tectonique 
verticale au point (x,y) et au temps t.. 

3.1.1. Les théories "globales" 

La première approche de la géomorphologie a été de considérer le système comme 
pouvant être décrit par une seule loi d'évolution reliant l'altitude au temps et à l'espace. 
Ces théories sont réductionnistes dans la mesure où elles réduisent l'ensemble des 
mécanismes potentiellement actifs à une seule loi d'évolution. Cette simplicité permet 
certaines prédictions, notamment sur la géométrie moyenne des profils. 

3 .1 .1 .1 . L évolution des pentes liées à des causes locales 

Dans son livre de 1961 ("Theoretical Geomorphology"), Scheidegger a résumé 
l'ensemble des lois simples, utilisables en géomorphologie, qui relient l'altitude et ses 
dérivées partielles (spatiales et temporelles) entre elles. Ces équations présupposent que 
les modifications morphologiques ne sont gouvernées que par les variations locales 
d'altitude. 

On peut alors représenter l'évolution des pentes par une équation différentielle du type : 



— = -V .q s +T(x,y,t) 
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où x est la distance suivant un axe horizontal, z l'altitude en fonction de x, et <P 
l'expression formelle de l'érosion. Scheidegger distingue trois types d'action érosive: 

(1) 4> = constante : le taux d'érosion est uniforme et indépendant de la pente; 

(2) <P = z : le taux d'érosion augmente avec l'altitude; par exemple lorsque la 
pluviométrie augmente avec l'altitude (pluies advectives); 

8z 

(3) 0 = 5— : le taux d'érosion est d'autant plus important que la pente est forte. 

ax 

d 2 z 

(4) <P = — - : le taux d'érosion est proportionnel au gradient de la pente, on retrouve une 

dx 2 

équation classique de diffusion où une perturbation de longueur d'onde X s'amortit avec 
un temps caractéristique te 

Scheidegger donne la solution de ces équations pour un versant, de pente initiale 
constante (Fig. II 13). Les solutions obtenues sont qualitativement raisonnables. 
L'expression de 4> peut être modifiée pour rendre compte de différentes formes 
plausibles pour les profils de versants. 

3.1.12. L'évolution des pentes prenant en compte l'aire drainée 

Les équations précédentes ne rendent pas compte d'effets non locaux qui peuvent être 
induits par la topographie amont. Parce que l'efficacité de la plupart des mécanismes 
d'érosion est dépendante du flux hydraulique, la loi d'érosion en un point est dépendante 
de l'ensemble des points en amont qui y drainent. Cette surface est désignée sous le terme 
générique d'aire drainée. 

Pour la plupart des mécanismes d'érosion, le flux de matière q s des systèmes 
géomorphologiques est usuellement exprimé en fonction de deux paramètres 
macroscopiques principaux, le flux hydraulique q et la pente S, avec une loi de transport 
du type : 

q s =F(q,S)=K q » S* 

où K est un coefficient. 

Citons quelques lois proposées dans la littérature : 
-le fluage ("creep") : n=0, m=l; 

-l'effet "splash" provoqué par l'impact des gouttes de pluie sur le sol et/ou le 
déplacement de blocs rocheux, n=0, m=2; 

-le transport de "sable" dans les chenaux alluvionnaires, n=\. 3-1.7, ra=1.3-2. 

La détermination de cette loi est principalement empirique. Pour certains processus tels 
que le ruissellement, il existe des justifications physiques basées principalement sur la 
contrainte cisaillante appliquée par la rivière sur son lit Le flux sédimentaire dépend de la 
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vitesse d'écoulement et de la géométrie de la rivière qui sont reliées à q et 5 via des lois 
statistiques de type équation de forme, ou équation de Manning ou de Chezy (donnant la 
vitesse du fluide). 

La résolution unidimensionnelle des équations de continuité et de transport (Smith et 
Bretherton, 1972) aboutit à un profil convexe lorsque l'exposant n est inférieur à 1 (par 
exemple dans le cas de la diffusion où n=0 et m=l), et un profil concave lorsque n>\. 
Ces auteurs montrent, de manière plus générale, que la condition de concavité implique 
l'inégalité suivante: 

dq Q 

Ces résultats analytiques confirment les travaux expérimentaux de nombreux auteurs 
sur les profils de pente (Culling, 1963; Kirkby, 1971, 1980, 1986). 

D'un point de vue purement formel, q le flux hydraulique est une variable externe qui 
possède sa propre équation de continuité : 

Vq = R(x,y,z,t) 

où R(x,y,z,f) est la pluie efficace (volume de pluie qui contribue au ruissellement) au 
point (x, y) et au temps t. 

Dans l'hypothèse où la pluie est uniformément distribuée (/?=c ste ), il est possible de 
remplacer la variable q par sa variable spatiale palliative: l'aire drainée, A. 

En utilisant la transformation q— »A, Howard (1994) illustre le rôle de l'exposant n sur 
la concavité de la pente, en prenant l'exemple de l'érosion d'un plateau dont la vitesse de 
surrection U est constante dans le temps et l'espace. Si l'équilibre érosion-surrection est 
atteint, on a : 

q s =U A= K A n S m 

soit, 

S=[U/K] 1/m A (1 - n)/m 

Si n>\, la pente diminue avec l'aire drainée, c'est à dire vers l'exutoire du BV; le 
système est donc concave. Sinon, la pente augmente avec l'aire drainée et la forme 
résultante est convexe. 

Pour conclure ce paragraphe, nous faisons trois remarques : 

(i) La conversion du flux hydraulique en aire drainée peut être une expression 

plus complexe que la simple loi: q~A. Des lois plus complexes, de type :q~À^ 
où y est un coefficient qui varie entre 0.6 et 1.0, ont été proposées pour rendre 
compte de configurations très différentes d'écoulement de surface. 
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(ii) Les équations de continuité et de transport sont vectorielles. Généralement, 
on suppose que les flux d'eau et de matières sont parallèles à la ligne de plus 
grande pente. Il est toutefois possible de définir une aire drainée qui 
corresponde pour des flux locaux divergents. 

(iii) Bien que systématiquement utilisée, l'assimilation flux hydraulique-aire 
drainée peut conduire à une solution inexacte 1) s'il existe une variabilité 
systématique de la pluie avec l'espace (en particulier une dépendance avec 
l'altitude) et 2) si les fluctuations temporelles de la pluie induisent des 
instabilités mécaniques. 

3.1.2. Les mécanismes couplés 

3.12.1. Les systèmes fluviatiles et la transition versant-réseau 

Il est évident que la formation de rivières modifie fondamentalement le processus 
d'érosion parce que le régime hydraulique est lui-même fortement modifié. En particulier, 
la vitesse hydraulique dans les rivières est de l'ordre de ~1 m-s" 1 , alors qu'elle est, en 
moyenne, plus faible d'un ou plusieurs ordres de grandeur sur les versants (-0.01-0.1 
m.s _1 dans les sols labourés). Le cisaillement appliqué sur les particules, proportionnel à 
la vitesse, est donc beaucoup plus intense sur les fonds de rivières que sur les versants. 

La modélisation des systèmes fluviatiles prend en compte deux propriétés principales : 

• une capacité érosive importante fonction du flux et de la pente 

• une capacité à transporter des charges solides dans la phase aqueuse 

Il est donc nécessaire de considérer au moins deux lois d'érosion distinctes pour 
différencier les processus lents et peu actifs agissant sur les versants et les processus 
rapides et efficaces agissant par l'intermédiaire des réseaux hydrographiques. 

Pour résoudre les équations d'érosion, il est donc nécessaire de définir une variable 
d'état qui détermine le processus actif en un lieu donné. Wilgoose et al. (1991a, b, c) 
introduisent une fonction "channel" qui possède sa propre loi d'évolution, fonction 
essentiellement de l'aire drainée locale. Ils reprennent des idées générales plus anciennes 
qui limitent l'existence de talwegs à une distance critique des lignes de crêtes (Horton, 
1945). Le lecteur trouvera dans les paragraphes suivants (Instabilités marginales en 
géomorphologie) une discussion plus approfondie de l'aire drainée critique et de sa 
signification physique. 

Cette dualité versant-réseau introduit une non-linéarité forte dans l'équation 
d'évolution du système mécanique liée à la migration spatiale de l'interface réseau- 
versant. 
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3.1 22. Le transport advectifdes rivières 

Les systèmes fluviatiles ont aussi la capacité de transporter des sédiments dans la 
phase aqueuse, en solution ou en suspension. Ce transport permet de déplacer rapidement 
de la matière sur de longues distances. 

La prise en compte de ce processus modifie considérablement la loi de continuité de la 
matière. Pour illustrer ce point, prenons l'exemple d'une érosion unidimensionnelle 
simplement proportionnelle à la pente. Sans transport par les rivières, l'équation de 
continuité est une équation de diffusion : 

<?L = D— 
dt djè 

Si, au contraire, le transport par les rivières est total, chaque élément érodé ne pourra 
pas se redéposer sur le point immédiatement voisin et l'équation de continuité deviendra: 

dz _ ^dz 
dt dx 

Dans le cas général, il faut considérer un transport partiel avec des échanges entre la 
surface topographique et la phase aqueuse. Pour cela, on fait intervenir une variable 
supplémentaire qui est la charge solide de la rivière, s. Le modèle le plus simple considère 
que l'augmentation de charge solide se fait proportionnellement à la pente locale, 5, et au 
débit hydraulique, q., tandis que la décharge (rivière— » surface) se fait 
proportionnellement à la charge solide. L'équation d'évolution de la charge solide 
devient 

= aqS-fis 

dt 

Si la pente 5 reste constante, cela conduit à une évolution de la charge s du type : 

Il existe donc (i) une charge limite a( j^ qui correspond à une stabilisation des 

échanges eau-sol, et une constante de temps II fi qui fixe la mise en charge de la rivière. 
On peut d'ailleurs remplacer ce temps par une distance caractéristique d=v/f5 où v désigne 
la vitesse des écoulements dans la rivière (Kooi et Beaumont, 1994). 

De même, l'équation de continuité de la matière doit prendre en compte les échanges 
sol-eau. On voit donc clairement que la prise en compte du transport advectif des rivières 
implique au moins une variable supplémentaire (la charge solide) et deux équations 
couplées. 
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3.1.2.3. l'altération chimique 

L'érosion d'un système géomorphologique est contrôlée par la capacité des sols à 
produire des matières transportables (l'altération, au sens large) et par leur capacité à 
transporter ce matériel mobilisable. On voit donc apparaître deux cas-limites : l'érosion 
contrôlée par l'altération ("Weathering-limited removal" ou "detachment-limited 
removal") et l'érosion contrôlée par le transport ("Transport-limited removal"). 

Peu de modèles géomorphologique prennent en compte l'altération comme facteur 
limitant de l'érosion. A notre connaissance, seuls Howard (1994a,b) ou Ahnert (1976, 
1977, 1987) l'ont fait. 

La méthode consiste à faire intervenir, en première approximation, une couche 
supplémentaire (la couche d'altérite aussi appelée régolithe) qui diminue par érosion via 
sa surface supérieure et augmente par altération mécanique ou chimique de l'encaissant 
via son interface inférieure. Ahnert considère que la vitesse de production des "altérites" 
dépend de l'épaisseur de la couverture régolithique, C, avec deux lois possibles : 

(i) pour l'altération mécanique, une loi exponentiellement décroissante, e"** 7 ; qui traduit 
la diminution des écarts thermiques avec la profondeur (effet de couverture); 

(ii) pour l'érosion chimique, une vitesse qui possède un maximum pour une épaisseur 
critique, C c . 

On a donc affaire à un mécanisme couplé où la vitesse d'altération dépend de 
l'épaisseur de la couche d'altérite, elle-même dépendante de l'érosion de surface donc de 
la loi de transport ! Dans le cas le plus simple, on peut faire intervenir une vitesse 
d'altération constante qui introduit un temps caractéristique, t a (inversement proportionnel 
au taux de production des altérite) dans le processus d'érosion. Une mise en situation de 
tels concepts dans un modèle d'érosion est présentée par Howard (1994), qui prend en 
compte le rapport d'érodabilité, F, entre le substratum et le régolithe. 

Une autre différence importante est faite entre le matériel mobilisé pour la première fois 
(substratum ou altérite) et les alluvions remobilisables dans le processus d'érosion. Ces 
derniers ont, a priori, des seuils d'érodabilité beaucoup plus faibles et un grano- 
classement déjà établi. Citons les études de Hairsine and Rose (1992a.) et Howard 
(1994) qui prennent en compte ce problème. 

3.2. Les instabilités en géomorphologie 

Il est évident que les processus géomorphologiques conduisent à la croissance 
d'instabilités mécaniques, géométriquement caractérisées par la formation de vallées plus 
ou moins profondes. L'étude de ces instabilités nous renseigne sur la physique de 
formation des reliefs. 
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Il existe deux façons diamétralement opposées d'aborder un tel problème 
(Loewenhertz, 1991). La première méthode est d'étudier la croissance d'instabilités 
marginales, c'est à dire l'évolution de petites perturbations dans un système proche de 
l'équilibre. Cette démarche a été entreprise par Smith et Bretherton (1972) et Loewenhertz 
(1991) entre autres. 

La seconde méthode consiste à établir de fortes perturbations dans les conditions aux 
limites du système (par exemple, la surrection brutale d'un plateau) et d'étudier la 
conséquence sur l'évolution du système. Cette méthode a été utilisée par Wilgoose et al. 
(1991), Stark (1994), entre autres. L'érosion associée est caractérisée par une phase de 
croissance forte des incisions. 

3.2. 1. La théorie des faibles perturbations 

L'évolution des systèmes géomorphologiques est fortement conditionnée par la 
croissance d'instabilités mécaniques, essentiellement liée au processus d'incision. 
Schématiquement, une incision, pratiquée localement, se traduit par une augmentation de 
l'aire drainée en ce point et/ou une augmentation de la pente locale. Si ces modifications 
d'ordres géométriques conduisent à une augmentation de l'incision par rapport au 
système environnant, on pourra parler d'instabilité mécanique. La condition dépend donc 
de la dépendance de la loi d'érosion F(A, S) aux fluctuations d'aire drainée A et/ou de 
pente 5. 

En termes mathématiques, il faut étudier la stabilité de l'équation de continuité de la 
matière pour des modifications faibles des paramètres d'entrée : z-»z+z', S->S+S', 
A-»A+A' où les perturbations d'altitude z', d'aire drainée A' et de pente S' sont petites 
devant les paramètres correspondants. Deux études ont particulièrement fait avancer ce 
problème : Smith and Bretherton (1972) et Loewenhertz (1991). 

Dans un système mono-dimensionnel, une incision ne crée pas de variation d'aire 
drainée. Seule l'augmentation de la pente locale modifiera la vitesse d'érosion. Pour tous 
les systèmes naturels, cela se traduit par une augmentation de la vitesse d'érosion et donc 
par diminution de la perturbation. Le système mono-dimensionnel est donc 
inconditionnellement stable. 

A deux dimensions, il existe des modifications d'aire drainée qui changent 
considérablement l'analyse précédente. L'apparition de flux transverses à la pente 
principale va avoir deux effets antagonistes : 

• favoriser l'érosion perpendiculaire à la pente principale (mécanisme stabilisateur), 

• concentrer ces flux transversaux dans les zones de bas-fonds dans une direction 
parallèle à la ligne de plus grande pente (mécanisme amplificateur d'instabilité). 



On peut montrer que, si le deuxième effet l'emporte sur le premier, l'amplification de 
perturbations transverses à la pente principale est possible. Dans tous les cas, les 
perturbations parallèles à la ligne de plus grande pente seront atténuées. 

Les analyses de Smith and Bretherton (1972) et Loewenhertz (1991) étudient la 
perturbation d'une solution stationnaire par une fonction sinusoïdale transverse de faible 
amplitude : z'{x,y,t)=W{x) exp (iky+ot) où x désigne la ligne de plus grande pente, y un 
axe perpendiculaire, t le temps, k le nombre d'onde de la perturbation (=2kIX où X est la 
longueur d'onde) et de facteur d'amplification. Si aest positif, la perturbation peut être 
amplifiée et donner naissance à une instabilité. Si G est négatif, la perturbation est amortie 
et le système est dit "stable". Notons que ces 2 études posent des conditions sur la 
solution de départ, avant perturbation, légèrement différentes, ce qui explique les 
différences dans l'analyse mathématique. Smith and Bretherton supposent que la solution 
de départ est de forme constante, c'est à dire que le taux d'érosion en chaque point est 
constant - en fait, leur solution est valable s'il n'existe pas, avant perturbation, de flux de 
matière transverses. Loewenhertz part au contraire d'une solution de départ réellement 
stationnaire où un soulèvement tectonique T , compense l'érosion en tout point. 

Les modifications d'aire drainée liées à la perturbation sont obtenues à partir de 
l'équation de continuité des flux hydrauliques. A partir des deux équations de continuité 
(matière et flux hydraulique), et en supposant que les flux se font le long de la ligne de 
plus grande pente, Smith and Bretherton proposent une condition nécessaire (mais pas 
suffisante) pour créer des instabilités: 

1s c dc ls 
q dq 

Malheureusement, comme le montre Loewenhertz (1991), cette condition n'est ni 
suffisante (une perturbation de longueur d'onde infinie est inconditionnellement stable 
puisqu'elle se ramène à une analyse unidimensionnelle), ni complètement nécessaire, bien 
que fortement recommandée. Smith and Bretherton montrent que la condition d'instabilité 
est équivalente à une concavité de la forme résultante stationnaire : 

& n 
— T <0 

La réciproque est aussi vraie. Ils montrent aussi que, si cette condition est vérifiée, les 
petites longueurs d'onde auront un facteur d'amplification a, proportionnel au carré du 
nombre d'onde, k 2 . Autrement dit, plus la longueur d'onde est faible, plus la perturbation 
sera amplifiée. 

Les conclusions de Loewenhertz (1991) sont qualitativement identiques. Toutefois, 
comme Loewenhertz part d'une solution stationnaire, la condition d'instabilité devient 
une condition de concavité de la surface de départ, avant perturbation. Elle pousse 
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l'analyse mathématique plus loin pour montrer que les instabilités peuvent quand-même 
croître au dessus du point d'inflexion, pour des points situés sur la partie convexe de la 
surface de départ. 

Ces analyses sont en désaccord, au moins sur un point, avec les observations 
géomorphologiques : ce ne sont pas les plus faibles longueurs d'onde qui sont les plus 
amplifiées. En effet, il semble qu'il existe une aire drainée critique en deçà de laquelle les 
incisions ne s'amplifient pas. Cette observation, souvent faite sur les systèmes naturels 
(voir Montgomery et Dietrich (1988, 1989)), a conduit Loewenhertz à modifier son 
analyse en introduisant une distance critique, y, en deçà de laquelle le flux ne se faisait 
plus suivant la ligne de plus grande pente mais de manière aléatoire tel un processus 
diffusif simple. Cette condition enlève le moteur de l'amplification des instabilités — la 
concentration des flux — pour des longueurs d'onde inférieures à y, et les longueurs 
d'onde préférentiellement activées ne sont plus les plus petites. La validité d'un tel 
modèle est liée aux mécanismes physiques qui justifient cette longueur y — effets 
hydrauliques inertiels ou gravitationnels, structure de la couverture végétale.... 
Raisonnablement, on peut penser que y ne dépasse pas la dizaine de mètre au grand 
maximum. Ce modèle ne peut donc pas être invoqué pour justifier des aires drainées 
critiques de l'ordre du km , comme celles observées sur le réseau hydrographique de la 
Vilaine (voir chapitre 1). 

Le concept d'échelle de coupure trouve déjà son origine dès les travaux de Horton 
(1945). Horton définit en effet une longueur critique d'écoulement laminaire de surface x c 
au-delà de laquelle la contrainte cisaillante t est supérieure à la contrainte cisaillante critique 
t c d'érosion. Des observations expérimentales (Dunne, 1980) montrent aussi qu'il existe 
une échelle critique en deçà de laquelle l'érosion ne produit pas d'incision mais se 
présente sous la forme de processus diffusifs, tels que l'effet splash (Fig. II 14). 

Ce concept d'aire drainée critique, qui définit la limite inférieure du réseau de drainage, 
a aussi été développé par Kirkby (1980, 1986), Tarborton (1992) d'un point de vue 
théorique, et par Montgomery et Dietrich (1988, 1989) d'un point de vue expérimental. 
La valeur de l'aire drainée critique est relative à un certain nombre de paramètres globaux 
ou locaux tels que: les propriétés locales du sol, les caractéristiques de la végétation, la 
pluviométrie ou encore la topographie locale. Les données expérimentales, obtenues sur 
trois bassins versants soumis à des pluviométries très différentes (260 mm/an, 760 mm/an 
et 1500 mm/an), révèlent une relation négative entre l'aire drainée critique a cr et la pente 
locale à l'amont de la source (Montgomery et Dietrich, 1988, 1989). Reprenant l'analyse 
de Tarboton et al. (1992), Montgomery et Foufoula-Georgiou (1993) construisent, à 
partir d'un MNT au pas de 10 m, la courbe générale du comportement de l'aire drainée en 
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b) 1/1500 sec 




c) 1/600 sec 




d) 1/300 sec 



d) 1/60 sec 



Fig 14: Effet de rejaillisement (splash). 

La dislocation des gouttes d'eau lorsqu'elles s'écrasent au sol 
provoque l'arrachement de particules de terres qui peuvent être 
projetées à une distance de quelques dizaines de c. 
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Figure H. 15 : Schematic illustration of relations between drainage 
area and local slope depicting hillslope/valley transition and channel 
initiation criteria. Note that the trend of the averaged data (wide 
shaded curve) for a catchment will vary between landscapes and 
watersheds. 

(d'après Montgomery et FouFoula-Georgiou, 1993) 



fonction de la pente locale à l'échelle du versant et des têtes de vallées (Fig II 15). 
L'utilisation d'un MNT à mailles triangulaires permet d'associer la partie croissante de la 
courbe aux zones divergentes (versant) et la partie décroissante aux zones convergentes 
(vallée). Ils concluent, contrairement à Tarboton, que le changement de pente ne 
correspond pas à un seuil de formation des cours d'eau pour un bassin versant donné, 
mais à la transition du versant à la vallée. L'existence de ces deux seuils pourrait 
s'expliquer par des constantes de temps très différentes entre les phénomènes 
hydrologiques et géomorphologiques. 

3.2.2. Les modèles déterministes de perturbations de fortes amplitudes 

32.2.1. Modèle de Willgoose et al. (1991) 

L'un des modèles continus de référence que l'on puisse citer dans ce domaine, est 
celui de Willgoose et al. (1990, 1991abc). Ce modèle utilise pour définir la dynamique du 
système, l'équation de continuité de la matière et une équation d'extension des réseaux: 



dt * y) Ps{l-n)\dx dxj [fol dy 2 ) 



dY_ 
dt 



0.0025 — o.irl 



(2) 



avec z l'altitude, t le temps, x et y les coordonnées horizontales, co le taux de surrection 
tectonique, r s la densité de sédiments, n la porosité, q s le flux de sédiments, D la 
constante de diffusion, Y le degré de formation du cours d'eau, dt le taux de croissance 
des cours d'eau, a la fonction d'initialisation des cours d'eau, at le seuil d'initialisation 
d'un cours d'eau. Les équations constitutives sont: 

a = p2q m2 S n2 

f(Y) = piO t Y = 0 (versant) 

A Y ) = Pi 7= 1 (cours d'eau) 

Qc=foA m > 
q = Q c /w 
w = p 4 Q c m 

où q est le débit d'eau, S la plus grande pente locale, b\ le taux de transport de 
sédiments, Ot le taux relatif de transport de sédiments dans le réseau, relatif à celui sur le 
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Séparation of trajectories with time ploned in log-log to 
illustrate the power-law growth of perturbations. 




Isometric view of channel network and hillslope with 
time: (a) t = 500, (b) t = 2000, (c) 1 = 6000, (d) t = 13,000. 
Times are nondimensional. 



Figures 11.16 : Simulation de croissance de réseau hydrographique à l'aide du modèle d'érosion de 
Willgoose: a), b) c), d); figures e) et f) étude de la croissance d'une perturbation en fonction du 
temps. Au- delà d'un certain stade de développement le réseau empêche toute nouvelle perturbation 
de croître, (d'après Willgoose et al. 1991a, Ijjasz-Vasquez et al., 1992) 



versant, Qc le débit local dans le réseau, w la largeur du cours d'eau, bi, mi et ni sont des 
constantes. 

Les termes de la partie droite de l'équation dynamique (1) représentent respectivement 
la surrection tectonique relative à l'altitude de l'exutoire, la continuité des transports de 
sédiments, et les processus de diffusion. L'originalité du modèle vient essentiellement de 
l'utilisation d'une fonction d'initialisation du réseau (2) empruntée à l'étude des 
croissances d'arborescences biologiques (Meinhardt, 1982). L'équation (2) possède deux 
attracteurs 7=0 (état initial du système, la surface est assimilée à un versant) et 7=1 (état 
stable irréversible, le pixel appartient au réseau). Lorsque la fonction d'initialisation a 
dépasse un seuil donné at , la valeur 7=0 devient instable et croit progressivement pour 
atteindre la valeur 7=1. Pour 0 < Y < 1, les équations 1 et 2 sont remplacées par: 



A Y ) = PiO t y < 0, la (versant) 



O, + {l-O t )(Y-0,la) 



1-0, la 

A Y ) = Pi 7>a (cours d'eau) 



0, la < 7 < a (état transitoire) 



où a est un paramètre supérieur à 1. Cette transformation est nécessaire si l'on veut garder 
la continuité du modèle. 

Le modèle de Willgoose est le seul modèle continu à notre connaissance qui différencie 
le réseau hydrographique du versant, d'un point de vue hydrodynamique. Les 
topographies numériques et les réseaux créés sont cohérents d'un point de vue qualitatif et 
quantitatif: lois de Horton et loi de Hack vérifiées. Willgoose montre, à l'aide de ce 
modèle, qu'une modification de la distribution aléatoire des altitudes initiales engendre, au 
cours de l'évolution du relief, des variations non négligeables de certaines propriétés 
géométriques et topologiques du réseau (lois de Horton et magnitude). Ce modèle montre 
également que plus les perturbations apparaissent tardivement au cours de l'évolution du 
système, moins elles sont amplifiées (Fig II 16). Pour un stade avancé du système, la 
perturbation régresse même pour ensuite disparaître. Les processus d'érosion semblent 
donc, d'après l'auteur de ce modèle, des processus déterministes hautement non-linéaires. 
De tels résultats concordent avec ceux des modèles analogiques (voir chapitre 2: "les 
modèles déterministes"). 

Si le modèle de Willgoose constitue jusqu'à présent l'un des efforts les plus 
conséquents en matière de modélisation déterministe, son nombre élevé de paramètres 
retire une grande partie de son intérêt 
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3.2.3. Les modèles stochas tiques de croissance 

Contrairement aux modèles déterministes basés sur des équations de continuité ou de 
transport, l'application des techniques de marcheurs ou de percolation dans la simulation 
des processus d'érosion-sédimentation ne présente pas de complexité excessive. Leur 
simplicité d'application favorisa l'éclosion d'un certain nombre de modèles, utilisant des 
règles de parcours, d'érosion et de sédimentation différentes (Takayasu et Inaoka 1992; 
Chase, 1992; Leheny et Nagel, 1993; Stark, 1994). Nous présenterons dans ce chapitre 
deux modèles que nous considérons parmi les plus achevés dans ce domaine. 

32.3.1. Modèle de Stark (1994) 

De manière à simuler l'érosion d'un plateau abrupt par effondrements successifs de 
son escarpement, Stark construit un modèle d'invasion percolation couplé avec de 
l'agrégation limitée par la diffusion. L'invasion percolation est en effet une technique 
relativement bien appropriée pour simuler le recul du front du plateau (escarpement) 
suivant le degré local (sur chaque pixel délimitant l'escarpement) de cohésion du matériau. 
La cohésion ou le niveau de dégradation sur chaque pixel ou chaque site de l'escarpement 
varie suivant l'exposition du site aux écoulements de surface, de subsurface et souterrains 
ou encore à l'érosion chimique. L'invasion percolation ne permettant pas de simuler les 
écoulements amont - aval, Stark propose d'introduire un modèle d'agrégation limitée par 
la diffusion pour simuler la dégradation des sites le long de l'escarpement au cours du 
temps. 

Les règles de parcours et d'érosion de ce modèle sont les suivantes: 

(1) sur chaque site du maillage représentant le système, on définit une valeur initiale 
aléatoire de résistance à l'érosion; 

(2) l'une des limites du système représente la position initiale de l'escarpement (le 
système est périodique sur les deux limites adjacentes au front de percolation; 

(3) chaque cycle d'érosion comprend 4 phases: 

(i) un marcheur aléatoire est parachuté de manière aléatoire sur un 
site, 

(ii) lorsque le marcheur atteint la limite du plateau en un pixel 
donné, la valeur de résistance à l'érosion de ce pixel diminue d'une 
quantité y, 

(iii) la résistance à l'érosion d'un site choisi aléatoirement sur le 
front de percolation (escarpement) est diminuée d'une quantité e, 

(iiii) le site appartenant au front de percolation et possédant la 
résistance à l'érosion la plus faible est érodé, à savoir retiré du 
plateau. 
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(4) lorsque le front de percolation atteint la limite opposée la simulation s'arrête. 
Suivant les valeurs de y et e , le modèle simule préférentiellement l'un des trois types 
d'érosion suivants: 

(i) e -> 0 et y-> 0 , le recul de l'escarpement est déterminé principalement par les 
conditions initiales sur l'ensemble du plateau (érosion d'invasion percolation), 

(ii) e » y » 0 , l'érosion à la base du plateau et à l'origine de l'évacuation des 
matériaux plus à l'aval, est prédominante et contrôle la progression du front de l'érosion; 
la dynamique des processus d'érosion est alors analogue à la progression d'un front de 
diffusion (érosion de diffusion); 

(iii) y» e » 0 , l'effet érosif des écoulements souterrains constitue la principale 
cause du recul de l'escarpement (agrégation limitée par la diffusion). 

L'érosion naturelle d'un plateau correspond à l'action combinée des trois types de 
processus précédents. Chacun de l'un de ces processus d'érosion engendre un front 
d'érosion (ligne irrégulière représentant dans un plan le bord du plateau) de dimension df 
différente. Cette dernière remarque s'applique également à la dimension d„ù n du chemin le 
plus court traversant la structure. Sur la base de cette différence, Stark propose un 
diagramme de phase tripolaire (Fig. II 17) dont chaque pôle symbolise un des 3 types 
d'érosion et correspond à un couple (df , d m i n ) donné. Le positionnement au sein de ce 
graphe de cas de figures réels, à partir de la mesure des dimensions df et d„à n , permettrait 
une estimation qualitative du processus d'érosion dominant. 

32.3.2. Modèle de Chase (1992) 

Conçu pour étudier l'évolution du relief sur des échelles de longueur et de temps 
importantes, le modèle de Chase est basé sur 3 règles simples simulant l'érosion, la 
sédimentation et la diffusion sur le parcours du marcheur. L'échelle du système est telle 
qu'il n'est pas nécessaire de distinguer les pixels appartenant aux réseaux 
hydrographiques de ceux des versants. Les règles de parcours et d'action du marcheur 
s'appliquent quelle que soit la localisation de ce dernier. Chaque marcheur représente un 
événement pluvieux suffisamment conséquent pour mobiliser des matériaux. 

Le marcheur se déplace suivant la plus grande pente et suit donc une règle de 
déplacement déterministe sans composante aléatoire. Les principales phases du parcours 
sont les suivantes: 

(1) parachuté de manière aléatoire sur la topographie, le marcheur provoque, lors de 
son impact avec la surface, une érosion amont-aval proportionnelle à la pente (coefficient 
d) entre le pixel d'impact et ses quatre plus proches voisins (pixels jouxtant l'un des côtés 
du pixel d'impact); cette première phase du parcours correspond, à l'échelle du système, 
aux processus de diffusion générés par l'altération de la roche au contact de l'eau, les 
glissements de pente, la formation de talus... ; 
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Invasion percolation 




Figure H. 17 : Eden g r °wth Laplacian fingering 

This is a kind of phase diagram intended to illustrate the relative importance of the 
model weathering and érosion processes. Invasion percolation, which reflects preexisting vari- 
ability in substrate strength and érosion rates, is unlikely to dominate érosion front morphology. 
Laplacian, or Poissonian, fingering arising from a DLA growth mechanism in the model reflects 
seepage weathering as well as seepage érosion and sapping and is thought to be a significant con- 
trol on plateau front morphology. Background weathering processes and mean substrate weakness 
are modeled through the Eden mechanism and are thought to be as important. Values of scaling 
exponents for the three end-member processes include d m i n , which is the fractal dimension of the 
shortest path across such structures, and which may affect the scaling of streams formed below 
the escarpment. Any mixed growth process may exhibit each scaling behavior over distinct length 
scale ranges, (d'après Stark, 1994) 



(2) une fois la phase d'impact terminée, le marcheur se déplace sur le pixel voisin (8 
voisins) le plus bas, et érode le pixel initial d'une quantité proportionnelle (coefficient é) à 
la différence de potentiel entre les deux pixels; le marcheur accumule ainsi une charge de 
sédiments, tant que sa capacité de transport n'est pas dépassée; 

(3) la capacité de transport du marcheur est également supposée proportionnelle à la 
pente; lorsque la pente décroît suffisamment pour engendrer une charge de sédiments 
dépassant la capacité de transport du marcheur, celui-ci dépose sur les deux derniers 
pixels de son parcours une quantité de sédiments telle qu'aucune inversion de pente n'est 
produite directement par la sédimentation; 

(4) le parcours du marcheur s'interrompt lorsque ce dernier atteint les limites du 
système ou un point bas de la topographie (les 8 voisins sont à des altitudes supérieures); 
le marcheur dépose alors toute la quantité de sédiments qu'il transporte. 

Une bonne compréhension du modèle nécessiterait plus de détails sur les règles 
d'action du marcheur, notamment au niveau des règles de diffusion et de sédimentation. 
Malheureusement dans son article de 1992, Chase ne détaille pas suffisamment son 
modèle. 

En jouant sur la valeur relative des coefficients d'érosion, de sédimentation et de 
diffusion, il est possible de mettre en évidence l'influence de ces trois processus sur 
l'évolution de la topographie. Afin de disposer d'un critère de comparaison, Chase 
quantifie les reliefs théoriques obtenus à l'aide du modèle en mesurant leur dimension 
fractale. Une topographie très irrégulière, où l'altitude en un point est indépendante de 
celles de ses voisins, présentera une dimension df tendant vers 3. A l'inverse, la 
dimension df d'un relief aplani tend vers 2. Le calcul de la dimension et l'étude qualitative 
des reliefs théoriques montre de manière très générale l'action des trois processus 
précédents. C'est ainsi que l'érosion advective accentue l'irrégularité des reliefs, la 
diffusion lisse les faibles longueurs d'ondes et la sédimentation atténue l'amplitude des 
grandes longueurs d'onde. 

Le modèle de Chase démontre qu'un modèle de marcheurs (3 paramètres) permet de 
construire des topographies réalistes en simulant de façon simple les trois modes de 
transport de sédiments que sont: la diffusion, les écoulements advectifs et la 
sédimentation. En ne distinguant pas le réseau hydrographique des versants, Chase 
suppose implicitement que les processus de diffusion et les processus d'advection 
agissent sur le même domaine spatial. Si le réseau hydrographique et le versant sont 
fortement couplés, ils ne peuvent cependant pas être confondus. Cette simplification 
introduit sans doute un biais non négligeable dans les résultats obtenus. 
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Les modèles d'érosion déterministes, basés sur la résolution des équations de 
continuité, nécessitent la définition de relations empiriques annexes et d'un nombre 
important d'hypothèses et de paramètres. Cette complexité contraste singulièrement avec 
la simplicité des modèles de marcheurs. Ces derniers, souples et performants, 
représentent une alternative de modèles d'érosion très intéressante à fort potentiel 
pédagogique. Les modèles de ce type, décrits dans la littérature, sont néanmoins d'une 
trop grande simplicité et ne distinguent pas, par des règles de marches, d'érosion et de 
sédimentation distinctes, les versants des réseaux hydrographiques. Il serait donc de 
notre point de vue intéressant d'améliorer les modèles de marcheurs. 

4. Un modèle couplé réseau-topographie 

Les analyses mécaniques précédentes montrent des approches des instabilités 
mécaniques très différentes. La première relie l'évolution du relief au développement 
d'instabilités marginales (faibles perturbations de la pente locale qui s'amplifient). La 
deuxième approche fait propager des perturbations à partir de fortes modifications de 
conditions aux limites du système. Aucun modèle actuel ne permet de relier ces deux 
approches du problème; sans doute parce que les simulations les plus réalistes des 
phénomènes d'érosion (Wilgoose et al. (1991), Stark (1992), Masek et Turcotte (1993), 
par exemple), considèrent implicitement que les réseaux ne peuvent se développer qu'à 
partir d'une des limites du système (par exemple, l'abaissement du niveau de la mer ou 
l'existence même de la mer aux bordures des systèmes continentaux). La croissance des 
réseaux représente alors une instabilité beaucoup trop forte pour permettre le 
développement d'autres instabilités marginales (instabilités de pente, par exemple). Ces 
modèles sont pourtant basés sur une ambiguité puisqu'ils définissent explicitement des 
règles de croissance liées aux conditions amonts (aire drainée, par exemple) alors 
qu'implicitement les conditions avales (existence d'un réseau qui draine vers l'extérieur 
du système) sont encore plus importantes. Bien entendu, les observations montrent que 
ces deux conditions sont remplies dans la plupart des systèmes naturels (sauf pour les 
systèmes purement endoréiques), mais on peut se demander si cette condition de 
connectivité-aval est une condition nécessaire du processus d'érosion ou une résultante 
(le problème de l'oeuf et de la poule, en quelque sorte !). 

Pour étudier l'influence respective des instabilités marginales et de la croissance à 
partir des conditions aux limites, nous avons construit un modèle probabiliste qui ne 
présuppose pas une condition de continuité des réseaux hydrographiques, tout en étant 
basé sur la dualité versant-réseau (liée aux conditions d'écoulement hydrographique) qui 
est, à notre avis, le moteur de la plupart des instabilités géomorphologiques. 
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La pluie représente l'une des sources principales djapport dans les systèmes 
géomorphologiques (l'autre source étant la tectonique). La variabilité spatiale et 
temporelle des précipitations se fait à des échelles de temps très courtes devant les 
vitesses de modifications du relief. Ce "bruit" de haute fréquence induit des fluctuations 
dans la réponse du système qui sont implicitement moyennées dans la plupart des 
modèles d'érosion. Nous avons aussi voulu modéliser ces fluctuations de haute 
fréquence et les conséquences qui en découlent sur l'évolution "moyenne" des variables 
d'état du système (altitude, aire drainée, flux hydraulique, nature des écoulements, ...). 

Concrètement, pour prendre en compte les fluctuations de haute fréquence des apports 
hydrauliques, nous avons défini la nature des écoulements hydrauliques (de versant ou de 
réseau) en utilisant un critère de flux instantané q plutôt qu'un critère portant sur sa 
variable géométrique palliative, l'aire drainée. En d'autres termes, nous voulons 
différentier le réseau de drainage (l'ensemble des points qui peuvent potentiellement 
devenir des réseaux hydrographiques) des points qui appartiennent réellement au réseau. 
La différence entre ces deux ensembles est essentiellement liée aux apports différentiels 
d'eau sur le système. Il est bien connu que les petites rivières s'assèchent en été... Avec 
ce modèle, nous définissons la probabilité qu'un point appartienne au réseau qui peut 
varier entre 0 (versant total) et 1 (réseau permanent). 

De plus, nous avons voulu construire un modèle où la topologie des connections ("qui 
draine dans qui !") n'est pas permanente mais peut varier, au moins à haute fréquence. 
La définition d'une aire drainée demande alors une moyennisation temporelle des flux. 
Avec les critères que nous utilisons, nous pouvons créer des flux divergents qui ne 
respectent plus la condition d'embranchement vers l'aval. 

Nous utiliserons ce modèle pour (i) étudier le rôle et l'importance des fluctuations de 
haute fréquence sur les variables d'état du système (état, flux, topographie, ...), (ii) 
caractériser les instabilités mécaniques 

4.1. Critère de formation des réseaux hydrographiques 

Les modèles de croissance de réseau utilisent dans leur grande majorité la notion d'aire 
drainée critique pour définir le passage du système versant au système hydrographique. 
La notion d'aire drainée est assimilée dans ce cas à la variable débit qui correspond 
beaucoup plus à la réalité physique du phénomène d'initiation des cours d'eau. 
L'utilisation d'un MNT, avec lequel seule une information topographique est représentée, 
ne permet pas dans l'absolu d'accéder à des variables dynamiques telles que les flux de 
matières ou d'eau. Ce n'est que par le biais d'une hypothèse d'apport d'eau uniforme sur 
l'ensemble du système, qu'il est possible d'interpréter les caractéristiques topographiques 
aux flux hydroliques. Cette hypothèse très restrictive, repose implicitement sur une 
échelle de temps suffisamment grande pour pouvoir moyenner les apports dans l'espace. 
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CRITERE: 

si (date courante - date dernier passage) < ou = temps critique 

-> formation d'un cours d'eau (création d'une maille de type réseau) 
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Fig 11.18: Critère d'initialisation d'un cours d'eau. 



Elle trouve sa justification dans l'existence d'une densité de rivières constante pour des 
zones de conditions climatiques uniformes. Le temps nécessaire à l'apparition d'une 
uniformisation des processus d'écoulement ou d'érosion varie selon le degré d'uniformité 
des conditions climatiques et pédologiques. Mais il n'a pas encore été démontré que 
l'antécédence de l'initialisation d'un réseau sur ses voisins était négligeable pour la 
répartition future des aires drainées. La chronologie des événements pluvieux locaux peut 
conditionner la prédominance de certains réseaux. Cet aspect n'est pas pris en compte 
dans les modèles qui confondent le débit avec l'aire drainée. 

Le modèle reprend le concept de précipitons introduit par Chase (1992) qui désigne 
des marcheurs dont les règles de déplacement sont fixées par la topographie locale, et qui 
ont une action de transport de matière. Les précipitons sont introduits aléatoirement sur la 
topographie. Le nombre de précipitons introduis définit le temps, ce qui revient à négliger 
le temps de parcours de chaque précipiton devant les fluctuations pluviométriques. 
L'unité de temps la plus petite, tu, (qui correspond à la plus haute fréquence de variation 
de la pluie) est l'écart de temps entre deux précipitons. Le flux hydraulique instantané est 
défini comme l'inverse de l'écart de temps entre deux passages successifs des 
précipitons, At, en un point de l'espace (que nous appelions "période de retour" dans la 
suite du texte). Cette condition revient à considérer que chaque précipiton transporte une 
quantité équivalent d'eau. Notons que cette valeur de flux est déjà une moyenne 
temporelle sur l'unité de temps tu* 

Nous ne définissons que deux types d'écoulement hydraulique possibles : un 
écoulement de rivière et un écoulement de versant Chacun de ces états hydrauliques va 
conditionner la loi locale de transport de matière du précipiton. Pour définir un critère 
d'état hydraulique qui prenne en compte les fluctuations de haute fréquence des apports 
pluvomiétriques, nous définissons un flux instantané critique, <7 m/ , en dessous duquel 
l'écoulement est de type versant et au-dessus duquel l'écoulement est de type réseau 
hydrographique. Avec les hypothèses formulées plus haut, ce critère revient à définir un 
temps critique, t cra , entre deux passages successifs deux précipitons : 

si At <î cr ih l'écoulement est de type réseau, 

si At >tcriu l'écoulement est de type versant 

Dans la pratique il nous suffira pour chaque maille m,y de garder en mémoire la date de 
passage du dernier précipiton (t re f), de la soustraire à la date du prochain passage d'un 
précipiton sur la maille (r re f - f c< > Mr ),pour en déduire la nature de la maille (si ( t re f - 
t cour )<t C rU alors rriij est réseau sinon m,y est versant), et d'actualiser la date mise en 
mémoire (Fig II 18). On prendra pour dater les passages des précipitons le nombre de 
précipitons lancés à partir du début de la simulation. 
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a) Grille des altitudes b) Grille des probabilités 

(Pi = (Ah|) 1/D ; D=coefficient de diffusion) 
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d) Déplacement du précipiton 

Fig 11.19 : Régie de déplacement sur le versant. 
Le déplacement s'éffectue suivant l'une des positions 
avales parmi les 4 directions possibles. On détermine 
pour chacune d'elles une probabilité Pi de déplacement. 
Le choix de la direction s'effectue en pointant de 
manière aléatoire l'une des position sur la longueur = EPi. 
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Fig II. 20 : Régie de déplacement dans le réseau. 
Le précipiton se déplace sur la position la plus basse 
par rapport à sa position courante parmi les quatres 
direction possibles 



On peut montrer que ce critère est équivalent à définir une aire drainée critique. En effet, 
la probabilité de recevoir un précipiton en une maille M quelconque de l'espace, pendant 
un temps t est : 

p=alA * t, 

où a est l'aire drainée au point et A l'aire totale du système (on prend l'unité de temps 
tu-=l). Statistiquement, le temps de passage moyen entre deux précipitons au point M 
correspond à l'inverse d'une aire drainée et le temps critique est donc représentatif d'une 
aire drainée critique avec une relation de type: 

a C rit=A* tultcrit- 

Contrairement au modèle de Wilgoose et al., les mailles n'acquièrent pas de manière 
irréversible le type "réseau". La fluctuation locale du nombre de passage de précipitons, 
due au lancé aléatoire, induit l'apparition temporaire de mailles de type réseau. Seuls les 
axes de drainage principaux gardent une certaine rémanence au cours de la simulation. 
Les variations les plus significatives de conversion versant->réseau ou réseau->versant 
sont situées principalement aux extrémités amont du réseau hydrographique, lorsque les 
bassins versants d'ordre 1 ne drainent pas suffisamment de précipitons pour maintenir 
une période de retour inférieure à la période critique. Une autre cause de tarissement des 
réseaux hydrographiques est l'apparition de zones planes situées à l'aval, où les taux de 
sédimentation sont supérieurs aux taux d'érosion, et où l'on observe des réorientations 
fréquentes d'écoulements avec l'apparition d'arborescences hydrographiques 
divergentes. Les arborescences divergentes engendrent une répartition diffuse des 
précipitons dans les différentes branches avales et par là même une augmentation de la 
période de passage. Au-delà d'un certain degré de divergence le réseau en terme de débit 
n'est plus viable et disparaît 

L'utilisation de la notion de débit à travers une fréquence de passage de précipitons, 
constitue l'une des originalités du modèle et permet sans nul doute de traduire avec plus 
de fidélité les phénomènes d'initialisation mais également de tarissement des cours d'eau. 
Ce dernier point est impossible à simuler en utilisant une variable palliative comme l'aire 
drainée critique. Nous pouvons donc à l'aide de ce critère distinguer, à chaque étape du 
parcours d'un précipiton, le versant du réseau, et appliquer les règles de parcours et 
d'action qui en résulte. 

4.2. Règles de marche 

Les processus de transport se divisent en deux classes: les processus diffusifs 
localisés principalement sur les versants et les processus advectifs dominants dans les 
cours d'eau. L'interprétation de ces deux types de transport en termes de déplacement de 
précipitons implique des règles différentes de marche pas à pas (Fig II 19 et H 20). Afin 
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différence de potentiel ( Ah^. ) 

Figure n.21 : Variation de la probabilité de déplacement en fonction de la différence de potentiel. 
Lorsque le coefficient de diffusion D augmente, la probabilité tend à se stabiliser autour de 1. 



de simplifier les calculs, nous utiliserons un déplacement de type D4, à savoir suivant 
l'une des 4 directions cardinales; le grands nombre de précipitons lancés aléatoirement 
compensant cette contrainte dans le choix des directions. Le déplacement est 
unidirectionnel, chaque précipiton étant considéré comme une entité indivisible. 

4.2.1. Sur le versant 

Les processus de diffusion purs sont caractérisés par un parcours aléatoire de 
particules dans un milieu donné. Leur modélisation se traduit donc dans l'absolu par une 
probabilité égale pour le précipiton de se mouvoir dans toutes les directions possibles. 

Cependant les processus d'érosion diffus comme les éboulis, les glissements de 
terrain ou l'effet splash sont soumis de façon non négligeable à l'action de la gravité. H 
nous semble donc important d'intégrer cette action dans le modèle. C'est pourquoi nous 
restreindrons le choix des directions possibles aux seules mailles ayant une altitude plus 
basse. 

On peut également accentuer l'effet de la gravité en introduisant la différence de 
potentiel Ah^ = (hi-hj) dans l'expression de la probabilité pij d'aller dans une direction j à 
partir de la position courante /. L'expression de pij en fonction de Ahij et d'un 
coefficient de diffusion D doit répondre à deux conditions: (i) elle doit être croissante 
lorsque Ahij augmente (plus la différence de potentiel Ahy est grande, plus les transferts 
de matières auront tendance à suivre cette direction), (ii) plus le coefficient de diffusion, 
D, est grand, moins le déplacement est influencé par la différence de potentiel Ahij. Nous 
avons choisi d'utiliser une fonction de probabilité en puissance de la différence d'altitude: 

Pij = -k Ah^ Y si Ahij>0 et /?,y=0 si A/*y<0; avec P tel que : pij=\ 



de manière à souligner le caractère diffusif des processus de transport sur les versants. 
La loi de puissance en 1/D atténue très fortement l'influence de la pente, lorsque D 
augmente (Fig II 21), et génère pour £>>10 des probabilités proche de 1 quelle que soit la 
différence de potentiel entre les points. Quand D-*0, on retrouve le critère de plus grande 
pente. 

4.2.2. Dans le réseau 

Les écoulements dans le réseau hydrographique sont supposés suivre la direction la 
plus grande pente et ne présenter aucune forme de diffusion. Leur modélisation est donc 
réduite à la simple détermination de la maille présentant le dénivelé le plus important par 
rapport à la position courante (Fig II 20). 
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Figure II 22a : La différence des taux d'érosion entre le versant et le réseau introduit une 
non-linéarité des processus d'érosion 



4.3. Règle d'érosion et de sédimentation 

La structure des reliefs témoigne de taux d'érosion de sédimentation sensiblement 
différents selon la nature physico-chimique de la surface et des processus d'érosion 
impliqués. Dans cette étude, on fera abstraction de la variété pédologique et géologique 
des substrats, pour se concentrer sur les conséquences de changements de processus 
d'érosion liées à la nature de l'écoulement 

Nous n'avons retenu que deux lois d'érosion possibles : 

-une loi de diffusion où la variation temporelle d'altitude est proportionnelle à la 
pente locale. 

-une loi de transport où il existe des échanges entre une phase aqueuse pouvant 
transporter des sédiments et le sol. 

4.3.1. Loi de diffusion 

Cette loi rend compte de déplacement locaux de la matière transportée. Traduit en 
terme d'action des précipitons sur la topographie, le transport local de sédiment 
correspond à la règle suivante: le précipiton érode une quantité de matière (CAhij) de la 
position courante /, proportionnelle à la pente locale (Ahij), pour ensuite déposer cette 
quantité de matière sur la position suivante j (Fig II 22a). La différence de processus 
entre réseau et versant est simulée en introduisant des coefficients d'érosion suivant la 
nature de la maille (réseau:CV ; versant:Cv). Dans l'approximation continue, cette loi de 
marcheur est équivalente à une équation différentielle du type: 

f = ^//(«V // A), 

où q désigne le flux hydraulique qui draine en ce point, et V// le gradient selon la ligne 
de plus grande pente (ou la moyenne pondérée des gradients selon les lois probabilités de 
déplacement définies plus haut). Notons que le flux drainé q dépend des règles de 
déplacement des précipitons. S'il n'y a pas de diffusion (D=0), ce flux est l'équivalent 
d'une aire drainée. Sinon, il faut définir une équation de continuité particulière du flux. 

Avec cette expression, la non-linéarité du problème mécanique n'est pas introduite par 
la loi de diffusion mais par la transition versant— >réseau qui correspond à un changement 
des paramètres C (C=C V pour le versant et C= C r pour le réseau). La loi de transport 
résultante, q s (q), est représentée sur la figure II 23. 

4.3.2. Loi de transport 

Les rivières possèdent une certaine capacité de transport qui dépend de la vitesse du 
fluide et de la granulométrie des sédiments mis en suspension (Howard, 1980). De 
manière générale une accentuation de la pente augmente le pouvoir érosif du cours d'eau 
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Etat final 



Erosion-sédimentation 

Fig ll.22bRègle d'érosion-sédimentation pour 
le versant et le réseau dans le schéma 1 et 
pour le versant dans le schéma 2 
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b) Erosion-sédimentation 

Fig II.23 : Règle d'érosion-sédimentation pour 

le réseau dans le schéma 2 

Cr: coefficient d'érosion dans le réseau, 

Ss: coefficient de mise en suspension 

Sd: coefficient de déposition, 

Stock: charge de sédiments dans le réseau 



et sa charge en sédiments. Inversement une diminution de la pente engendre un dépôt de 
sédiments et une réduction de la charge de matière en suspension. Le cours d'eau se 
comporte comme un réservoir de sédiments qui possède sa propre loi d'évolution. Le 
modèle d'érosion-sédimentation consistera donc à compléter le schéma précédent en 
ajoutant une variable supplémentaire : le stockage de sédiments dans la phase aqueuse, s. 
Bien entendu, cette loi de transport ne s'applique qu'aux mailles de type réseau. Les 
règles qui suivent ne concernent donc pas les mailles de type versant, pour lesquelles les 
règles d'érosion-sédimentation (diffusion) restent inchangées. 

Nous utilisons une loi de transport simple où l'augmentation de s est proportionnelle à 
la pente locale, et où une proportion s/Ld de la charge solide est redéposée, Ld étant 
l'analogue d'une longueur caractéristique (voir le paragraphe sur le mécanismes couplés 
au début du chapitre 2). 

Lorsqu'un précipiton se déplace d'une position i à une position j, il érode la position i 
d'une quantité Ss Ahij, qui est stockée dans la phase aqueuse. Simultanément, le 
précipiton redépose au point j , une proportion s/Ld de la charge en sédiments (Fig 
II 22b). Cette redéposition est toutefois limitée pour empêcher une inversion des pentes 
lors de la phase de sédimentation. Si le précipiton quitte le réseau pour aller sur des 
mailles "versant", il continue à déposer des sédiments suivant les règles précédentes, sur 
autant de mailles que son stock le lui permet. Précisons que lorsque le précipiton tombe 
dans un trou, il libère tout le stock emmagasiné avant de disparaître. Ce cas de figure peut 
générer une inversion de la pente locale si le stock dépasse le volume du trou. Cette 
modification du relief reste néanmoins très localisée et n'engendre pas de perturbation 
sensible dans l'évolution du relief. 

Dans la limite continue, ces règles de transport sont analogues aux équations décrites 
dans le paragraphe 3.1.2.2. (le transport advectif des rivières). Ld représente une 
distance caractéristique du transport 

Si besoin, cette équation de transport peut se combiner simplement avec une équation 
de diffusion. 

4.4. Résultat du modèle de précipitons 

L'utilisation de ce modèle très schématique vise à mettre en évidence l'influence des 
principaux types de processus d'érosion sur le développement des réseaux 
hydrographiques et du relief. Nous distinguons dans le modèle deux modes de transport 
de matière: l'un diffus et l'autre advectif, et deux groupes de processus suivant leur 
localisation sur le versant ou dans le réseau. Le modèle interprète sous la forme d'un 
coefficient, le taux d'érosion ou la capacité de transport caractéristique de chacun de ces 
sous-ensembles. En jouant sur les valeurs relatives des coefficients, il ainsi est possible 
d'accentuer ou d'atténuer l'érosion différentielle engendrée par le changement de régime 
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hydraulique. Un rapport — - élevé traduit un fort pouvoir d'incision du ruisseau, et une 

c v 

valeur élevée le transport sur de longues distances de quantité de matière érodée 
localement avec un taux de transport proportionnel à S s Ld . 

L'échelle de temps des processus est implicite suivant la valeur des coefficients des 

paramètres C r et C v . Entre deux passages de précipitons l'altitude en un point varie en 

fonction de la valeur de ces deux paramètres. Considérons deux simulations, une 

première réalisation avec des valeurs de coefficients d'érosion de référence, et une 

C 

deuxième réalisation avec des valeurs de Cr et Cv. supérieures, mais avec un rapport 

C y 

identique. Les variations d'altitudes à chaque pas de temps seront supérieures dans le cas 
de la deuxième simulation. Le rapport constant des coefficients d'érosion, induit 
implicitement que les processus sont identiques dans les deux réalisations. Seule une 
érosion sur une plus longue durée permet alors d'expliquer des altitudes finales 
inférieures. Le temps écoulé entre deux lancés de précipitons n'est donc pas le même 
pour les deux simulations. Les processus d'érosion ne sont définis dans le modèle que 
par la valeur relative des coefficients d'érosion. Des valeurs différentes mais un rapport 
constant ne changent en rien le cas de figure traité. 

4.4.1. Evolution des réseaux hydrographiques et des reliefs. 

Nous présentons dans ce paragraphe deux exemples de simulations afin d'illustrer les 
principales phases d'évolution de réseaux de drainage. Ces deux exemples correspondent 
à des échelles de longueur et de temps différentes de plusieurs ordres de grandeur. La 
première simulation reproduit l'installation d'un réseau de rigoles sur un versant de pente 
constante et de rugosité initiales égale à 1% de l'altitude maximale du relief. Le deuxième 
exemple simule l'érosion progressive d'un plateau d'élévation moyenne constante 
(z(f=0)=l). 

La limite aval du système Cc=0) est maintenue à une altitude nulle; à l'amont (x=l), le 
flux de matière est nul. Des conditions périodiques sont imposées sur les limites latérales 
du système. La principale différence entre les deux simulations vient de la perturbation 
initiale imposée par la condition aval. Dans le cas de la première simulation, il n'existe 
pas de différence d'altitude entre la limite aval du système et la base du versant. Dans le 
cas du plateau, il existe initialement une forte différence d'alitude entre l'extrémité du 
plateau et la limite aval. 

4.4.1.1. Erosion d'un versant de pente constante 

Nous posons pour cette simulation les hypothèses suivantes: 

(1) les processus de transport dominants sur le versant et dans le réseau sont 
diffusifs (pas de stockage de matière dans le réseau: 5s = 0); 
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(2) les écoulements sur le versant ne sont pas suffisamment importants (charge 
hydraulique des filets d'eau faible) pour induire une érosion différentielle 
importante entre les écoulements concentrés (rigoles) et les autres processus 
érosifs dominants sur le versant (-=r=2). L'effet de dispersion est caractérisé 

C y 

par un coefficient de diffusion D=l. 

Le temps critique /^est égal à 20 de manière à laisser se développer un réseau dense 
de rigoles. 

Deux types d'information, l'une topographique (blocs diagrammes) et l'autre 
hydrographique (réseau hydrographique dans un plan), sont présentés. A tout instant, un 
pixel peut avoir soit un état Réseau, soit un état Versant. Nous avons calculé, pour 
chaque pixel, la probabilité qu'il soit à l'état Réseau, c'est à dire le nombre d'occurrence 
où il est dans cet état normalisé au nombre d'occurrence où le pixel a été considéré. Les 
figures montrent plusieurs stades d'évolution correspondant chacun à N précipitons 
lancés (N=15 000 * numéro de l'étape d'érosion). 

Les premiers stades d'évolution (figure II 24) présentent une organisation rapide du 
réseau de drainage. Les précipitons ne suivent pas au stade 1 de parcours bien défini. 
Mais dès le stade 2 des chemins préférentiels arborescents apparaissent. Le réseau de 
drainage ainsi initié s'organise par la suite (Fig n 24 et n 25, stades 3 à 25) par une 
succession de captures, en un réseau stable ne présentant plus de variations de structures 
significatives (Fig II 26). La croissance du réseau de drainage se fait par coalescence de 
structures initialement isolées et non pas par croissance progressive de structures 
arborescente à partir de la condition aux limites avale. Dans cet exemple les instabilités à 
l'origine de l'incision du versant ne sont donc pas induites par des modifications de 
conditions aux limites du système. L'organisation du réseau correspond plutôt au 
développement d'instabilités marginales. Il est intéressant de noter que le réseau acquière 
une relative stabilité dès les premiers stades. Le relief au contraire suit une évolution plus 
linéaire avec une amplification progressive des structures topographiques (Fig II 25 et II 
26). 

4.4.1.2. Erosion d'un plateau 

Les hypothèses utilisées sont les suivantes: 

(1) le transport de matière dans les cours d'eau est purement advectif (Cr=0); 

(2) le rapport des coefficients d'érosion entre la diffusion de versant et le 

S 

transport , est égal à 4; 

C v 

(3) la longueur caractéristique du transport, Lj, est égal à 20 pixels, soit 1/6 de 
la largeur du système. 
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Figure II 24: Erosion d'un versant de pente constante; premiers stades d'évolution du 
réseau. 




Figure II 25: Erosion d'un versant de pente constante; premiers stades d'évolution du 
relief. 




Figure II 26: Erosion d'un versant de pente constante; phase de stabilisation de la 
structure des réseaux. 




Figure II 27: Erosion d'un versant de pente constante; stades ultérieur d'évolution du 
refief. 



Le temps critique t cr it et le coefficient de diffusion D sont inchangés (t cr i t =20, D=l). 
Notons cependant que le nombre N de précipitons lancés entre chaque étape présentée est 
10 fois supérieur au précédent (N=150 000 * numéro de l'étape d'érosion). 

Cet exemple montre essentiellement 3 phases d'évolution (Fig II 28 et II 29): 

(1) une croissance aval-amont du réseau hydrographique avec la formation de 
lacs dans les zones planes (stade 10) 

(2) un développement maximal des réseaux hydrographiques avec une 
connexion de tous les points au niveau de base (disparition des lacs) (stade 20); 

(3) la formation de zones de dépôt (vallées alluviales) engendrant des 
arborescences hydrographiques divergentes vers l'aval (les cours d'eau 
deviennent moins localisés). 

Contrairement à l'exemple précédent le réseau ne possède pas de structure stable et 
évolue fortement avec le relief. La variabilité des réseaux théoriques obtenus et la 
possibilité de créer des divergences Amont/Aval mettent en valeur le potentiel du modèle à 
reproduire une évolution réaliste des réseaux hydrographiques. Les critères de validation 
courant tels que les lois de Horton ou de Hack ne peuvent être utilisées pour de telles 
configurations. Ces dernières n'ont en effet de sens que dans le cas de réseaux 
arborescents convergents d'Amont en Aval. Ceci retire une grande part de l'intérêt que 
pouvaient représenter ces lois, quoique d'utilisation déjà critiquable (voir Chapitre I § 7 et 
Chapitre II § 2.2.2.2). Nous nous contenterons donc d'une analyse qualitative des 
résultats. 

La croissance des réseaux simulée à l'aide de ce modèle n'obéit pas à un mode unique 
de progression. L'extension des réseaux s'effectue à la fois par une progression des têtes 
de biefs vers l'amont, et par la capture de structures non raccordées au niveau de base. 
Cette simulation montre aussi le développement de structures de drainage non connectées 
à la limite aval du système. Indépendantes des conditions limites du système, elles ne 
sont contrôlées que par les configurations locales de la surface. Ces structures finissent 
par être capturées par le réseau principal . Ce schéma de progression est sensiblement 
différent, et à notre sens plus réaliste, de celui d'une coissance pas à pas du réseau à 
partir de la limite aval du système, qui était proposé jusqu'à présent dans les modèles. 

Les deux exemples de réalisation démontrent que l'application simple de règles 
d'érosion et de sédimentation suffisent à reproduire un certain réalisme. A défaut d'une 
analyse quantitative, l'évolution des réseaux théoriques obtenus a permis de faire ressortir 
l'une des principales originalités du modèle: la propagation des réseaux ne peut se réduire 
à une croissance irréversible vers la condition aux limites aval. Les phénomènes de 
capture et les instabilités de pente jouent un rôle prépondérant sur l'évolution du relief 
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Figure II 28: Erosion d'un plateau; stades d'évolution du réseau. 




Figure II 29: Erosion d'un plateau; stades d'évolution du relief. 



Notre critère de réseau, en intégrant la variabilité locale de la fréquence des pluies et en 
supprimant les conditions de branchement aval, procure une certaine liberté d'évolution 
aux réseaux vis à vis des conditions topographiques. Dans les modèles utilisant une aire 
drainée critique, la géométrie des réseaux est entièrement contrôlée par la topographie. 
Une modification du réseau est forcément liée à une modification topographique. Le 
critère de débit critique tel qu'il est conçu, à savoir avec une forte composante aléatoire 
due à la pluie, introduit un degré de liberté supplémentaire au système. 

Ce modèle démontre de plus qu'il n'est pas nécessaire d'introduire une hypothèse 
implicite de connexion comme un critère d'aire drainée amont pour engendrer des réseaux 
arborescents. 

4.4.2. Analyse quantitative des instabilités développées par les processus 
d'érosion 

L'analyse qui suit vise à définir l'influence de paramètre sur l'évolution du profil 

moyen du système et sur la croissance d'instabilités. Les paramètres étudiés sont: le degré 

de diffusion, l'érosion différentielle entre le versant et le réseau, ou encore la capacité de 

transport des cours d'eau. Plusieurs simulations ont été effectuées en utilisant des 

r 

combinaisons des paramètres (D, t cr it , jr-, Ld, Ss) différentes. Pour chacune d'entre 

elles une analyse fréquentielle par Transformée de Fourier permet de dégager l'évolution 
des principales longueurs d'ondes et par là même des instabilités qui se développent 

Le profil initial choisi pour toutes les simulations est un profil sinuisoïdal, 
convexe/concave, avec pour condition limite amont (x=l) un flux nul, et pour condition 
limite avale (x=0) le maintient d'une altitude nulle. Des conditions périodiques sont 
imposées sur les bords latéraux. La condition limite-aval implique une évacuation totale 
de la matière arrivant en x=0. La rugosité de la topographie initiale est égale à 1% de 
l'élévation maximale. 

4.42.1. Erosion de versant, influence de paramètre de diffusion D 

L'influence du paramètre D sur le développement des instabilités est testée en ne 

r 

considérant que le transport sur les versants (77-= L Ld,=Ss =0). Les processus d'érosion 
sont donc supposés linéaires et purement diffusifs. 
Trois cas sont étudiés: 

- D=0 correspondant à des déplacements suivant la plus grande pente locale 
quelle que soit la position du précipiton; 

- D=l où la probabilité d'aller dans une direction est proportionnelle à la pente; 

- D=10 pour lequel les directions de déplacement sur des points avals sont 
pratiquement équiprobables. 
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Les topographies résultats (Fig H 30, II 31, II 32) présentent un profil moyen évoluant 
d'un profil convexe/concave à un profil concave. On observe en effet une migration du 
point d'inflexion vers l'amont, ce dernier disparaissant par la suite progressivement. 
Cette tendance vers un profil concave trouve son origine dans l'application de 
déplacements diffusifs orientés vers l'Aval. Les règles de déplacement sur le versant 
lorsque D^O ne correspondent pas à une véritable diffusion isotrope. Ce cas de figure ne 
peut donc pas être assimilé au cas théorique d'un profil convexe/concave soumis à des 
processus diffusifs isotropes. Rappelons que la partie convexe d'un tel profil soumis à 
des processus diffusifs devrait en théorie être stable. 

Quelle que soit la valeur de D, les instabilités développées sur la partie convexe du 
profil sont atténuées dans la partie concave. On peut en déduire que le modèle d'une 
érosion diffusive linéaire ne permet pas aux instabilités de se développer dans un profil 
concave. Il est donc nécessaire, pour développer des instabilités dans un profil concave, 
d'introduire une non linéarité dans la modélisation des processus érosifs. 

Suivant la valeur du coefficient de diffusion, l'orientation des structures 
topographiques est plus ou moins contrainte par le pendage moyen du profil. Lorsque les 
déplacements s'effectuent suivant la plus grande pente (D=0), l'orientation des talwegs 
suit la pente moyenne du système. L'introduction d'une composante aléatoire dans 
l'orientation des déplacements génère une dispersion des directions de drainage 
notamment dans la partie concave du système. 

4.422. Les instabilités liées à la formation des réseaux 

r 

L'influence du rapport des taux d'érosion R--^~ peut être mis en évidence en 

imposant un coefficient de diffusion fixe (D=l), et une capacité de transport nulle dans le 
réseau (Ld,=Ss =0). Le temps critique de formation des cours d'eau t cr it est égale à 20 de 
manière à laisser la possibilité au réseau de se développer. Nous avons simulé l'évolution 
du profil convexe/concave pour trois valeurs de R (R=l, R=4, R=10). Une analyse de 
fourrier sur 63 stades d'évolution permet de suivre dans le temps la valeur et l'amplitude 
des longueurs d'ondes dominantes pour chacun des cas de figure. 

L'étude qualitative des topographies résultats (Fig II 33 et II 34) montre qu'une 
augmentation du rapport R engendre une accentuation du relief et l'installation plus rapide 
d'un profil concave. La différence de potentielle entre les points situés dans le réseau et 
ceux situés sur les lignes de crêtes augmente d'autant plus que le pouvoir d'incision du 
réseau est important. Nous observons alors une augmentation des pentes et de ce fait une 
accentuation de la rugosité du relief. L'analyse des longueurs d'ondes dominantes sur 
une ligne du profil (ligne 100, ligne où l'amplitude des longueurs d'onde est maximale) 
révèle (Fig H 35 et H 36): 
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Figure II 30: Influence des processus de diffusion sur l'évolution d'un versant 
convexe/concave. Cas D=0: les écoulement suivent la plus grande pente. 




Figure II 31: Influence des processus de diffusion sur l'évolution d'un versant 
convexe/concave. Cas D=l: écoulements seirri-détenriirrites vers l'Aval. 




Figure II 32: Influence des processus de diffusion sur l'évolution d'un versant 
convexe/concave. Cas Z)=10: les écoulements s'effectuent de manière aléatoire vers 
l'Aval. 



Figure II 33: Influence de la non linéarité des taux d'érosion entre le versant et le réseau 
sur l'évolution d'un profil convexe/concave (pas de transport advectif de matière dans le 

réseau). R=^ =4 



Figure II 34: Influence de la non linéarité des taux d'érosion entre le versant et le réseau 
sur l'évolution d'un profil convexe/concave (pas de transport advectif de matière dans le 

réseau). R=&- =10 
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Figure II 35: Variation en fonction du temps de l'amplitude de longueurs d'onde 
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dominantes et du rapport des taux d'érosion R=^f- > (pas de transport advectif de matière 
dans le réseau). 
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Figure II 36: Variation en fonction du temps de l'amplitude de longueurs d'onde 
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dominantes et du rapport des taux d'érosion R=^r- , (pas de transport advectif de matière 
dans le réseau). 



- des variations d'amplitude plus marquées pour des valeurs importantes de R, 

- la nécessité d'imposer une érosion différentielle pour développer des 
longueurs d'onde significatives, 

- une variation de la longueur d'onde dominante en fonction du temps. 

La décroissance, de l'amplitude des longueurs d'ondes dominantes pour i?=10, 
s'explique par une diminution sensible des variations d'altitudes, qui favorise une 
diffusion plus importante et par conséquent un lissage des longueurs d'ondes. On peut 
distinguer ainsi 3 phases dévolution: 

- une phase d'installation d'une longueur d'onde caractéristique, pendant 
laquelle la longueur dominante décroît pour atteindre une valeur stable; 

- une phase de stabilité; 

- une phase d'aplanissement de la topographie, au cours de laquelle on observe, 
la disparition progressive du relief due à une dominance des processus 
diffusifs, et l'apparition d'une pénéplaine. 

L'existence de longueurs d'ondes bien définies dans l'évolution du relief démontre 
que la non linéarité des taux d'érosion entre le versant et le réseau est une condition 
suffisante pour permettre une croissance des instabilités. Notons que l'amplitude des 
longueurs d'onde est maximale à l'amont de la transition convexe/concave Nous 
observons ainsi encore une fois l'indépendance des phénomènes de croissance de 
perturbations vis à vis des conditions aux limites du système dans le cas d'un profil 
convexe/concave. 

4.4.2.3. Influence du temps critique tcrit déformation de cours d'eau 

Quatre cas de figures ont été étudiés de manière à mettre en évidence l'influence du 

débit critique de formation des réseaux hydrographiques. Les figures II 37, II 38 et II 39 

correspondent à des temps critiques de 10, 20, 40. 

Lorsque le temps critique augmente (ou l'aire drainée critique diminue), les longueurs 
d'onde de la topographie tendent à diminuer. Cette tendance s'explique par 
l'augmentation de la densité du réseau avec ce paramètre t C rit • Pour d es aires drainées 
critiques faibles, on observe donc bien une forte densité de réseau, et inversement. 

Pour t C rit =5 nous n'observons pas d'instabilité sur le profil convexe/concave. Dans ce 
cas, la pérennité des réseaux n'est pas suffisante pour induire des instabilités. L'érosion 
est dans ce cas dominée par la diffusion de matiière sur les versants. Avec les mêmes 
conditions mais en utilisant une géométrie initiale de type plateau (Fig II 40), nous avons 
pu développer des instabilités en incision. Ceci montre que des perturbations fortes aux 
limites du système peuvent être suffisantes pour développer des instabilités mécaniques. 
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Figure II 37: Influence du débit critique nécessaire à la formation d'un cours d'eau sur 
l'évolution d'un profil convexe/concave, (pas de transport advectif de matière dans le 

réseau) t cr it =10 



Figure II 38: Influence du débit critique nécessaire à la formation d'un cours d'eau sur 
l'évolution d'un profil convexe/concave, (pas de transport advectif de matière dans le 

réseau): t cr it =2 0 



Figure II 39: Influence du débit critique nécessaire à la formation d'un cours d'eau sur 
l'évolution d'un profil convexe/concave, (pas de transport advectif de matière dans le 

réseau): t cr it =4 0 




Figure II 40: Influence du débit critique nécessaire à la formation d'un cours d'eau sur . 
l'évolution d'un plateau, (pas de transport advectif de matière dans le réseau): t C rit=5 



4.42.4. Les instabilités liées au transport advectifpar les rivières 

Pour décrire l'influence des phénomènes de transport advectifs de matière dans les 
réseaux sur la croissance des instabilités, deux cas de figures ont été envisagés: l'un 
présentant une forte différence des taux d'érosion entre le versant et le réseau 
hydrographique (Fig II 41 et II 42), et l'autre se limitant à au cas d'une érosion égale 
dans le réseau et le versant mais avec une capacité de transport dans le réseau importante 
(Fig II 43 et II 44). A la différence des simulations sans transport de matière, ces deux 
exemples montrent des structures très marquées liées à la non linéarité induite par les 
phénomènes de transport. La cinétique de l'évacuation de la matière dans le réseau est 
plus importante que dans le cas simple de transports uniquement locaux (diffusion sur le 
versant). L'effet lié à l'érosion différentielle entre versant et réseau est amplifiée par le 
transport de sédiments dans la phase aqueuse. On note cependant que les phénomènes de 
transport ne suffisent pas au développement d'un réseau arborescent (Fig II 41), mais 
permettent néanmoins l'installation d'un réseau bien défini sans structures divergentes 
(Fig Et 43). La non linéarité induite par les phénomènes de transport n'est donc pas une 
condition d'initiation mais d'amplification. 
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Figure II 41: Influence du transport advectif de matière dans le réseau sur l'évolution 
d'un profil convexe/concave (pas d'érosion diférentielle entre le versant et le réseau): 
configuration du réseau. 




Figure II 42: Influence du transport advectif de matière dans le réseau sur l'évolution 
d'un profil convexe/concave (pas d'érosion diférentielle entre le versant et le réseau). 




Figure n 43: Influence du transport advectif de matière dans le réseau sur l'évolution 
d'un profil convexe/concave avec une non linéarité des taux d'érosion entre le versant et 
le réseau: configuration 



Figure II 44: Influence du transport advectif de matière dans le réseau sur l'évolution 
d un profil convexe/concave avec une non linéarité des taux d'érosion entre le versant et 
le réseau. 



5. conclusion 

Une approche statique de la structure des réseaux hydrographiques, telle que celle 
développée au chapitre 1, ne permet pas de remonter aux principes physiques 
d'organisation, n est donc essentiel de percevoir le caractère dynamique de ces systèmes 
très évolutifs, à travers une modélisation de leurs phases de croissance. Une telle 
approche nécessite une définition des principaux phénomènes d'érosion à l'origine de la 
transformation des reliefs et de l'installation des cours d'eau. L'analyse des processus 
d'érosion montrent que les mécanismes d'altération d'une part et les mécanismes de 
transport d'autre part ont des cinétiques différentes. Une représentation simplifiée conduit 
à un schéma non linéaire dans lequel on distingue deux sous systèmes: (i) un réseau, où 
les taux d'érosion sont importants et les transports rapides et advectifs; (ii) un système 
complémentaire que nous appelons "versant", sur lequel les taux d'érosion sont moindres 
et les transports locaux et diffus. 

L'utilisation de ce schéma dans un modèle de croissance des réseaux hydrographiques 
permet d'appréhender l'influence des principaux paramètres physiques sur le 
développement des instabilités à l'origine du relief et du réseau hydrographique. Nous 
avons, dans cette perspective, développé un modèle basé sur le concept d'entités 
élémentaires ("précipitons") représentant un volume d'eau donné. Les précipitons se 
déplacent sur la topographie et la modifient par une action érosive. L'emploi d'entités 
mobiles permet d'étudier l'effet des fluctuations de haute fréquence des apports d'eau sur 
le développement du système. 

Nous avons pu ainsi mettre en évidence que les instabilités dynamiques de la 
topographie pouvaient être induites soit par l'amplification d'instabilités marginales dues 
à des perturbations locales des pentes, soit par des perturbations liées aux conditions aux 
limites. Le relief est donc la résultante de la croissance simultanée de ces deux types 
d'instabilités. On peut néanmoins, en jouant sur les variables d'état du système (altitude, 
aire drainée ...), favoriser l'un ou l'autre. On montre aussi qu'une loi d'érosion non 
linéaire par rapport aux flux hydrauliques est une condition nécessaire au développement 
d'instabilités. Le couplage avec un transport dans les rivières est un facteur 
d'amplification. En revanche, les processus de diffusion tendent à inhiber le 
développement de ces instabilités. 



150 



CHAPITRE 3: GÉNÈSE DES RÉSEAUX HYDROGRAPHIQUES, 
APPROCHE EXPÉRIMENTALE 



1. Introduction 

L'étude quantitative de la structure des réseaux hydrographiques à l'échelle du Massif 
Armoricain, et les résultats encourageants des modèles d'érosion, ont mis en avant la 
nécessité de comprendre les processus de formation des cours d'eau. La précision avec 
laquelle l'extension du réseau hydrographique est déterminée tient en effet une place 
importante dans la compréhension du développement des systèmes hydrographiques et 
topographiques. Des caractéristiques géométriques et une localisation géographique 
précise des biefs d'ordre 1 permettraient sans nul doute de mieux définir l'ensemble des 
caractéristiques géométriques du réseau. Les biefs extérieurs constituent en effet la moitié 
de la population de biefs d'une arborescence binaire. Une meilleure estimation de leur 
longueur conduirait dès lors à une définition plus précise de l'échelle de coupure en deçà 
de laquelle le réseau hydrographique ne suit plus une organisation aléatoire. Mais c'est 
avant tout un suivi dans le temps du déplacement des points sources qui apportera et 
autorisera une appréciation plus juste des critères de formation. Une vision statique du 
réseau hydrographique ne peut correspondre à la réalité de systèmes évolutifs tels que les 
cours d'eau. La compréhension de l'organisation des réseaux et une estimation des temps 
de transferts caractéristiques passent donc par l'étude de la dynamique des sources dans 
le paysage. 

Cependant dans le cas de régions humides, bien des sources ne sont pas ponctuelles 
mais diffuses et se présentent sous la forme de surfaces saturées en eau. Leur extension 
varie en fonction des épisodes de pluie, on parle ainsi de zone à surface variable. Ainsi 
l'étude de la variation des surfaces saturées connexes au cours d'eau est nécessaire dans 
l'analyse de l'évolution du réseau hydrographique selon les conditions d'humidité. 

Localisation des points sources 

Très peu d'études expérimentales ont apporté des indications sur la localisation des 
sources en fonction des caractéristiques physiques ou géométriques du relief. Les travaux 
de Montgomery et Dietrich (1988, 1989) mettent cependant en évidence une corrélation 
entre l'aire drainée amont et la pente locale mesurée à des points sources appartenant à 
trois bassins versants différents (Fig. III. 1). Une mesure de ces deux variables 
topographiques, à l'amont et à l'aval des points sources, se traduit par la distinction de 
deux populations de points de part et d'autre de la courbe définissant la transition de l'état 
versant à l'état ruisseau. L'extension des réseaux hydrographiques serait donc contrôlée 
par un seuil défini par des variables topographiques. Comme le fait remarquer Kirkby, 
(1988) l'analyse de Montgomery et Dietrich, porte sur un échantillon de sources 
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Figure m.l : Mesure de l'aire drainée en fonction de la pente locale le long des lignes de talweg. Les 
graphes A B et C représentent le rélevé effectué sur 3 bassins versant différents. Les points marqués 
d'un triangle représentent une mesure effectuée à l'aval de la source, les ronds pleins au niveau de la 
source, les points ajourés à l'amont de la source. Les données des points sources définissent un seuil 
de transition topographique entre versant au réseau (D). (d'après Montgomery et Dietrich, 1992). 




Figure III. 2 : Les extrémités amont des biefs 
d'ordre 1 peuvent se présentés sous la forme 
d'un effondrement local (a) ou sous la forme 
de zones saturées diffuses. Dans le deuxième 
cas de figure, il est impossible de définir 
géographiquement la source par un point, 
(d'après Kirkby, 1988). 



localisées à des incisions très nettes dans le sol et par là-même ne recouvre pas l'ensemble 
des configurations possibles de formation de cours d'eau. Bien des ruisseaux se 
présentent à leur extrémité amont sous la forme de zones diffuses saturées en eau ou de 
segments discontinus de petits biefs (Fig. III.2). La détermination d'un point 
géographique précis de formation du cours d'eau s'avère dès lors pour le moins difficile. 
L'image d'un réseau fîlaire bien délimité par ses points sources est en réalité une vision 
très simpliste et ne peut pas rendre compte de tous les cas de figures rencontrés dans les 
régions humides telles que la Bretagne. De plus cette vision puremenent spatiale ne tient 
pas compte de l'influence de l'état hydrique du bassin versant et du régime 
pluviométrique dans l'extension du réseau; deux variables qui font implicitement 
référence à une évolution temporelle du milieu. 

Les processus hydrologiques à l'échelle du versant impliquent des échelles de temps 
beaucoup plus faibles que celles jusqu'à présent considérées dans ce travail. L'extension 
vers l'amont ou la régression vers l'aval des biefs sources, pendant une pluie ou un cycle 
saisonnier, peuvent être en effet significative. La dépendance de la position géographique 
des sources vis-à-vis d'événements pluviométriques ou de l'état hydrique initilal est 
d'autant plus marquée que le milieu présente un pouvoir de diffusion (végétation, 
rugosité du sol, perméabilité du substrat) et une capacité de stockage des volumes d'eau 
faibles. Sous des climats arides des réseaux importants se créent au cours d'une pluie ou 
d'une période pluvieuse pour disparaître lors de périodes sèches. Sous des climats 
tempérés et humides, où les réseaux sont moins éphémères, le réseau rémanent est utilisé 
comme structure de référence dans les modèles. Mais signalons que même en Bretagne 
des biefs d'ordre 1 peuvent présenter des variations de longueur de l'ordre du kilomètre 
au cours d'un cycle saisonnier. La représentativité du réseau rémanent n'est donc que très 
relative et insuffisante pour traiter les problèmes liés à la réponse transitoire d'un bassin 
versant à une précipitation. En toute logique une étude expérimentale rigoureuse 
nécessiterait un suivi de l'évolution du réseau à l'échelle du versant pendant et entre 
chaque pluie. Une telle étude demanderait dans l'absolu une photographie de l'ensemble 
du bassin versant d'ordre 1, avec un pas de temps inférieur d'un ordre de grandeur à la 
durée de la précipitation. Les mesures de la position des sources, présentées dans la 
littérature (Montgomery et Dietrich, 1988, 1989), sont généralement plus représentatives 
du débit de base que du débit de pointe. Plus généralement la mesure et l'interprétation de 
la position des sources dans la topographie, impliquent donc d'être liées à des échelles de 
temps hydrologiques (pluie, crue, saison, année), pour lesquelles on détermine une 
localisation moyenne pour la durée d'observation. 
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Figure III. 3a: Séparation des hydrogrammes de crue par une méthode 
isotopique. L'eau préexistante à l'averse peut représenter jusqu'à 90% des 
débits instantanés de crue, (d'après Merot, 1981). 



1.2. Extension de la notion de réseau hydrographique 

La dichotomie réseau/versant est fondée sur la distinction de deux comportements 
dynamiques, l'un caractérisé par des mouvements d'eau et de matière rapides et 
unidirectionnels (m/s ou cm/s), et l'autre par des transferts lents et diffus (10" 6 m/s). Si 
l'on s'attache à distinguer les compartiments géographiques d'un bassin versant selon 
leur réponse vis-à-vis d'un apport d'eau, la différenciation versant-réseau n'est pas 
appropriée aux régions humides où le rôle des zones saturées est important 

Le concept de zone contributive fut introduit pour la première fois par Cappus (1960) 
et repris par la suite par les auteurs anglais Hewlett et Troendle (1975). Il désigne sous un 
terme générique les surfaces du bassin versant participant au transfert rapide de volumes 
d'eau dans le réseau. La vision trop simpliste d'un ruissellement par refus à l'infiltration, 
homogène sur l'ensemble du bassin versant (Horton, 1933), n'intègre pas l'observation 
d'écoulements préférentiels sur les zones saturées en eau. D'autre part, l'application des 
méthodes isotopiques de séparation des hydrogrammes, révèle un apport important d'eau 
antécédente à la pluie pendant la montée de la crue (Fig. m.3a) (Crouzet et al., 1970; 
Pilgrim et Huff, 1978; Sklash et Farvolden, 1979; Merot et al., 1981; Pearce et al., 1986; 
Sklash et al., 1986). Ce phénomène s'explique par la mise en charge des horizons de sol 
saturés en eau ce qui provoque, par un effet de chasse, un flux d'eau préexistante dans le 
drain aval à l'origine de l'essentiel du pic de crue (Fig. m.3b). On considère ainsi que le 
ruisellement est non pas généralisé sur tout le bassin versant avec un taux variable suivant 
la perméabilité, mais limité aux zones contributives avec un taux de ruissellement de 100 
% (Fig. m.4). Dans un tel schéma, la réponse rapide du bassin ne se localise plus dans le 
seul réseau hydrographique mais également au niveau des surfaces saturées. Si l'on 
différencie les unités fonctionnelles d'un bassin versant suivant le temps de transfert qui 
les caractérise, les zones saturées ne sont pas dissociables du réseau hydrographique. 
L'extension du réseau est fortement dépendante de l'évolution parallèle de la surface 
occupée par les sols saturés. La compréhension des processus de formation des cours 
d'eau passe donc par une étude du couplage entre le réseau et les zones contributives. 

La dynamique des zones saturées, pour des conditions initiales d'humidité données, 
est essentiellement contrôlée par la topographie et le type de sol (Beven, 1978; Beven et 
Kirkby, 1979; O'Loughlin, 1981). Cependant deux critères ont orienté rapidement 
l'étude des variations des zones saturées vers une approche purement topographique. 
D'une part la disponibilité et la qualité de l'information topographique sont très 
supérieures à celles de la répartition spatiale des sols, et d'autre part, bien qu'il n'existe 
pas encore de modèle explicite, la typologie des sols est en partie liée à des critères 
topographiques (Beven, 1978; Mérot et Ezzahar, 1991; Widiatmaka, 1994). On peut 
donc espérer à travers l'analyse du relief local retrouver la répartition spatiale des 
occurences de saturation du sol. Le rapport al fi où a représente l'aire drainée par unité de 
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Figure in.3b : Contribution des écoulements souterrains à la génération des écoulements 
rapides de crue: la mise en charge des profils de sol engendre par un effet de chasse 
(piston flow) une alimentation rapide du cours d'eau en eau préexistante à l'averse. 
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Figure IIL4 : Deux schémas possibles de fonctionnement d'un bassin versant, (source- 
Mérot, 1994). 



contour de courbe de niveau, et fi la pente locale, constitue la structure de base des 
indices proposés dans la littérature (Beven et Kirkby, 1979; O'Loughlin, 1986). Il 
correspond à l'intuition que le temps de résidence des volumes d'eau augmente lorsque la 
pente décroît, et que l'occurence de saturation est proportionnelle à l'aire drainée amont 
Plus la valeur de l'indice associé à un élément du bassin versant est élevée, plus cet 
élément tendra à être saturé en eau. Un tel indice ne tient compte que de la géométrie du 
sous bassin situé à l'amont du point considéré, en faisant abstraction du niveau de base 
de la nappe. L'étroite relation entre les cours d'eau et les zones saturées qui leur sont liées 
laissent néanmoins supposer une influence avale non négligeable des axes de drainage 
sur l'extension de la saturation des sols, n nous semble donc intéressant d'étudier la 
position relative des surfaces saturées par rapport aux axes de drainage, et par là-même de 
tenter de mieux définir l'extension des points sources en terme d'augmentation des 
surfaces contributives sur le bassin versant. 

2. Étude expérimentale de la distribution spatio-temporelle de 
l'humidité de surface 

2.1. Les approches choisies 

L'humidité des sols dépend de nombreux paramètres tels que les caractéristiques 
physique du sol, la topographie, la végétation, les séquences climatiques locales... Elle 
présente ainsi une variabilité très forte tant sur le plan spatial que temporel, et nécessite un 
grand nombre d'échantillons afin d'obtenir des estimations fiables de distributions 
spatiotemporelles. La répartition spatiale des zones contributives ou plus généralement de 
l'humidité des sols est par ailleurs difficile à modéliser dans l'espace. C'est pourquoi peu 
de travaux traitent de l'évolution de la répartition spatiale de l'humidité des sols sur des 
échelles de l'ordre d'un bassin versant d'ordre 1 (Gascuel-Odoux, 1984). La difficulté à 
mesurer l'humidité in situ ou à la modéliser, stimula les recherches vers l'utilisation de 
techniques de télédétection, qui offrent la possibilité d'acquérir une information sur de 
grandes surfaces avec une précision de l'ordre de la dizaine de mètres. 

L'étude expérimentale présentée suit deux approches: 

(1) le suivi, à l'aide d'images radar enregistrées tous les 3 jours, de la variation, à des 
échelles de temps proche de l'averse, des surfaces saturées en eau sur toute une région 
(100 km x 100 km); l'objectif final étant de pouvoir obtenir une statistique représentative 
de la variation des zones de sources, en observant un grand nombre d'ordre 1; 

(2) des campagnes mensuelles d'analyses d'échantillons de sol sur un bassin versant 
d'ordre 1, de manière à suivre l'évolution de la distribution spatiale de l'humidité de 
surface au cours d'une saison. 
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Précisons que l'approche (1) nécessite une première phase de validation de 
l'information radar sur un petit bassin versant. Nous avons donc, pour mener à bien les 
deux approches, choisi un bassin versant test d'ordre 1, situé sur la trace au sol suivie par 
le satellite. Le site expérimental est un sous-bassin versant du Coët-Dan (1,3 km 2 ) 
affluent de l'Evel dans le Morbihan. Il présente un relief peu prononcé sur du schiste 
briovérien. La pluviométrie annuelle moyenne est de 711 mm/an avec des apports 
journaliers rarement supérieurs à 20 mm/ j. Les sols se développent dans une couverture 
limoneuse sur du schiste briovérien. L'origine des limons vient d'une part de l'altération 
des substrats géologiques en place et d'autre part de remaniement éolien. Le bassin 
versant du Coët-Dan est le siège d'une culture intensive (92% de la surface du bassin sont 
exploités) et d'une production de porcs élevée (4 fois la moyenne bretonne). 

2.2. Application de la technique radar 

2.1.1. Mesure de l'humidité par technique radar 

Contrairement aux fréquences dans les domaines de l'infrarouge et du visible, les 
micro-ondes (1 - 300 GHz) possèdent des caractéristiques essentielles dans la perspective 
d'une mesure de l'état hydrique d'un bassin versant. Elles sont peu sensibles à l'état de la 
couverture nuageuse; elles ont une capacité plus ou moins grande de pénétration sous le 
couvert végétal et dans les premiers centimètres du sol; et elles présentent un 
comportement variable suivant les propriétés diélectriques du sol, que l'on sait dépendre 
de la teneur en eau. La possibilité de mesurer l'humidité moyenne des sols à l'aide de 
micro-ondes actives fut démontrée théoriquement et expérimentalement dés les années 70 
(Ulaby et Batlivala, 1976; Ulaby et al., 1978, 1983; Bernard et al., 1982, 1986; Bruckler 
et al., 1988). Il ressort de ces travaux que le signal rétrodiffusé de l'émission radar, est 
fortement corrélé avec rhumidité des 10 premiers cm de sol. Notons également une forte 
dépendance de la valeur du signal avec l'état de surface du sol (rugosité, orientation des 
structures principales) et de la couverture végétale qui absorbe et diffuse tout ou partie les 
ondes radar. Les paramètres instrumentaux peuvent être néanmoins choisis de manière à 
minimiser les perturbations induites par l'état de la surface ou du couvert végétal: une 
fréquence proche de 5 GHz (bande C), un angle d'incidence compris entre 10° et 20°, et 
une polarisation HH. En utilisant une configuration du radar adéquate, il est possible 
d'étalonner le signal rétrodiffusé en terme d'humidité volumique (Ulaby et al., 1978; 
Bernard et al., 1982, 1986; Bernard et Vidal-Madjar, 1983) (Fig. E1.5). Les étalonnages 
proposés dans la littérature ne sont valables cependant que pour des valeurs humidité 
inférieures à la capacité au champ. Lorsque les sols atteignent des taux d'humectation 
supérieurs, la réponse en terme de signal rétrodiffusé n'est plus linéaire. L'apparition 
d'une couche d'eau libre en surface induit un effet spéculaire et une chute brutale du 
signal (Fig. m.6). Les courbes expérimentales d'étalonnage précédentes deviennent alors 
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Figure m.5 : Calibration du radar Erasme (signal rétro diffusé en fonction de l'humidité des sols 
entre 0-10 cm). Remarque: les humidités mesurées ne dépasse pas la capacité au champ, (d'après 
Jarry, 1987). 
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Figure 111*6 : Direction et intensité du 
signal rétro diffusé en fonction de l'état 
de la surface. Pour une surface lisse (par 
exemple recouverte d'eau) le signal émis 
vers le radar (direction des ondes 
incidentes) est nul; suivant la rugosité de 
la surface la diffraction de l'onde 
incidente est plus ou moins isotrope, 
(source: Trevett, 1986) 
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caduques, du fait de l'impossibilité de pouvoir distinguer le signal diffusé par un sol très 
sec, de celui obtenu par la réflexion sur un sol saturé. Ce comportement singulier permet 
toutefois de déterminer la localisation des sols présentant une humidité supérieure à la 
capacité au champ dans un contexte climatique humide (Brun et al., 1990). H est donc 
permis d'espérer délimiter l'extension des zones saturées à partir d'images radar sur de 
grandes surfaces. 

2.1.2. Les images radar du satellite ERS 1 utilisées - Résultats 

Le radar SAR (Synthetic Aperture Radar) embarqué sur le satellite ERS 1 (mis en 
orbite en Mai 1991) enregistre le signal rétrodiffusé sur de 100 km de côté avec un pixel 
de 12,5 m x 12,5 m. Le programme d'observation de la Terre du satellite ERS 1 
comporte essentiellement deux phases de cycle de répétitivité orbitale: une phase de 
périodicité à 3 jours pour l'observation de sites particuliers, une autre avec une périodicité 
de 35 jours pendant laquelle l'ensemble de la surface terrestre est couvert. La situation 
privilégiée du bassin versant du Coët-Dan (12 km 2 , Morbihan), situé sur l'une des traces 
au sol suivi par le satellite pendant la phase à 3 jours, a permis l'acquisition d'images 
radar sur ce site tous les 3 jours de Janvier à Mars 1992, et tous les 35 jours pendant 1 
an. Nous pouvions donc tenter de suivre l'extension des zones saturées sur plusieurs 
bassins versants d'ordre 1 pendant toute la période hivernale. La configuration du SAR 
embarqué sur ERS 1 est la suivante: 5,3 GHz, VV, 23°. Un contrôle au sol de l'humidité 
est nécessaire pour valider l'étalonnage. Une série de campagnes de mesures de 
rhumidité de surface, et un suivi de l'extension des zones saturées sur un des sous 
bassin versant d'ordre 1 du bassin du Coët-Dan ont été effectués. 

L'analyse de la variation du signal à l'intérieur de chaque parcelle confirme la bonne 
sensibilité du signal à l'humidité de surface. Elle montre néanmoins une très forte 
dépendance de la valeur du signal en fonction de la rugosité de la surface des sols et du 
couvert végétal. Cette dépendance est très marquée dans le cas de cultures saisonnières où 
l'évolution de la couverture végétale introduit des variations significatives de l'intensité 
du signal rétrodiffusé (Fig. IU.7). Il est donc difficile de traduire en terme d'humidité 
l'information acquise à 35 jours d'intervalle (intervalle pendant lequel on observe 
d'importantes modifications du couvert végétal et de la rugosité du sol, surtout dans le 
cas de régions agricoles telles que la Bretagne). Pour les images obtenues avec une 
périodicité de 3 jours, il est possible d'atténuer l'effet dû à une modification de la 
géométrie de la surface de diffusion, en étudiant seulement la différence de l'intensité du 
signal entre deux images successives. La non linéarité du signal n'est observée que pour 
quelques cas de figures, à savoir: une récolte sur la parcelle, un travail du sol, le gel ou 
l'extension des zones saturées. Ce traitement permet de montrer pour certaines parcelles, 
une corrélation entre une chute de l'intensité et la saturation des sols (Mérot et al., 1994; 
Appendice). Mais le bruit important généré par l'état de la surface nécessite un traitement 
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Figure ni .7 : Suivi du signal radar par type de culture. L'intensité du signal rétrodiffusé 
suit des variations caractéristiques du couvert végétal. 



du signal radar plus conséquent pour acquérir une information utilisable sur toute 
l'image. Une telle analyse d'image demande une technicité en traitement du signal radar et 
un suivi simultané des caractéristiques de la surface du sol ou du couvert végétal, qui sort 
du cadre de cette thèse. C'est pourquoi il n'a pas été possible d'utiliser la télédétection 
radar pour suivre l'extension des zones contributives, à l'échelle d'un bassin versant tel 
que le Coèt-Dan puis d'une région plus large. 

2.2. Approche topographique de la distribution spatio-temporelle 
de l'humidité de surface 

Les campagnes d'échantillonnage, effectuées en premier lieu dans la perspective d'une 
validation terrain des images radar, ont été analysées dans un second temps afin d'extraire 
la distribution spatiale de l'humidité en fonction de la topographie, au cours d'un cycle 
hydrologique. 
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ABSTRACT 



Time and space distribution of water content at soil surface has been investigated on a small 
farmland catchment (1.3 km 2 ) from four field surveys corresponding to différent moisture status. For 
each survey, about 400 samples were collected at soil surface and at 5 cm in depth along ten axis 
parallel to the greatest slope. The relationships between thèse measurements and the topography are 
analysed. The structure of the data is well explained by a topographie index referring to the 
downslope conditions and defined as the élévation différence between the sample point and the 
stream point corresponding to the outlet of the water pathway derived from Digital Elévation Model 
(DEM). This index can be considered as an hydraulic head, at least in saturated conditions. A 
threshold on this index allows two domains in the catchment to be distinguished; the water content is 
nearly constant and slowly variable in the upper domain while it strongly increases and is highly 
variable in the bottom domain. A decreasing exponential function of water content versus this index 
fits well with the data averaged by index class. The physical significance of this relationship for the 
catchment processes, implications for sampling stratégies of remote sensing validation are 
discussed. 



WORDS Water content Space/time distribution Radar calibration Topographie index 
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INTRODUCTION 



A major challenge of hydrology or bioclimatology is the assessment of time and space distribution of 
soil surface water content at the catchment scale, that would greatly improve hydrological modelling, 
particularly for distributed and physically based models using the contributive area concept 
(Anderson and Kneale, 1980; ; Band, 1986; Beven and Kirkby, 1979; Beven and Wood, 1983). 
O'Loughlin (1986) then Barling et al. (1994) proposed a prédictive model of spatial distribution of 
saturated areas and soil water content based on this concept: a wetness index is estimated by the 
model as only the DEM and a value of transmissivity are required a priori. Expérimental data or 
estimâtes based on expérimental data of soil surface water content would be useful to validate 
modelling assumptions or to better constrain initialization and calibration of an hydrological 
modelling. 

As a matter of fact, expérimental studies are scarse while the processes are numerous and 
occur at différent time and space levels. At the catchment level, soil surface water content dépends 
on the location in hillslope, due to the water table. Locally, it dépends on the physical characteristics 
of soil surface and layers, and on surface conditions such as roughness or overlying végétation. 
Many informations are so required to study experimentally the spatial distribution of soil surface 
conditions, i.e. a detailed topographie analysis, a delineation and a characterization of parcels and 
soils considered as relevant units. 

To make up for the lack of expérimental studies the active microwave remote sensing (Bernard 
et al., 1982; Brun et al., 1990; Ulaby and Batlivala, 1976; ) particularly by using the data of the 
spaceborne synthetic aperture imaging radar (Mérot et al., 1993) may be considered as a real mean 
to estimate the spatial distribution of soil surface moisture. However the ability of the backscattering 
signal to provide it really dépends on expérimental calibrate phases. Even in this way, it is primary 
necessary to investigate how the soil surface moisture is distributed on the catchment in order to 
choose available sampling stratégies of calibration and to get unbiaised spatial estimâtes. 

The purpose of this paper is to analyse the influence of the topography on the spatial distribution 
of soil surface water content, on a small catchment of 1 .3 km 2 from four field surveys corresponding 
to différent moisture status on the catchment. The relevance of a model and the conséquences on 
calibration phase of a remote sensing approach are discussed. 

EXPERIMENTAL SITE 

The study catchment is an area of 1 .3 km 2 that corresponds to a subcatchment of the Coët-Dan 
catchment (French Brittany) (Fig. 1). The orientation of the main stream is East-West. The south- 
facing slope présents an élévation range of about 20 mètres; it is an open field landscape of fodder 
and cereal crops. The north-facing slope présents a lower élévation range of about 10 mètres; it is 
mainly covered by grassland. The bedrock is comprised of primary brioverian shales. Soils are loamy 



brown low leached soils, sometimes stoneless. The fairly homogeneous texture is the following: loam 
70%, clay 17%, sand 13%. They are classified in three groups according to their drainage régime: no 
waterlogged, waterlogged in depth, waterlogged from the surface. The areas of thèse three soil types 
are similar ; the most waterlogged areas are located along the channels. There is a high contrast of 
saturated hydraulic conductivity from one to two order of magnitude from well drained slope soils to 
the waterlogged bottom soils (Widiatmaka and Curmi, 1994). This catchment présents nether roads 
nor hedges. 

MATERIAL AND METHOD 

Fourfield surveys were performed on 18 February 1992, and 24 October 1992, 27 January 1993 and 
3 September 1993. The cumulative rainfall for the previous month is respectively 52, 38, 105 and 3 
mm. Thèse surveys ranges from the wettest to the driest conditions observed during 1992 and 1993: 
in january 1993, the survey took place in very wet conditions whereas in september 1993 it took 
place in very dry conditions. 

The sampling design is comprised of ten axis parallel to the greatest slope to study specifically 
topographie effects (Fig. 1): the greatest slope is assumed to contrai the soil moisture distribution 
along the hilislope. The axis are distributed throughout the catchment to include equally the three soil 
groups and ail types of land use. The sampling lag on each axis is of 25 mètres; about 20 samples by 
depth were collected on each axis. The soil was sampled at two depths: 0 to 5 cm, and 5 to 10 cm. 
500 samples were collected for each survey except for the last one: only one axis out of two was 
sampled in september 1993. Each sampling point is referred by its geographical coordinates 
determined using a 1/10000 scaled topographie map. This map is also used to compute the DEM 30 
mètres grid size that allows to détermine the élévation of each sampling point. 

Gravimétrie water content was determined on each sampie of each survey by weighting it before 
and after48 hours oven drying. A few samples obtained from saturated areas présent free water and 
gravimétrie water content are over-estimated. Thèse outlayers influence further caiculations; 
therefore, thèse samples were set to a maximum water content of 80% but not eliminated to 
incorporate information on the location of the saturated areas. Bulk density was determined by 
cylinder method for one sampling point out of three at 5-10 cm during one survey. 

In a first stage, spatial distribution of soils, land uses are considered in order to analyse the 
spatial distribution of thèse measurements. The three soil units delineated according to their drainage 
régime with a 1/10000 scaled soil map are distinguished. Four land uses are distinguished: bare soils 
or soils with low végétation cover; soil surface with more than 50% of végétation cover; plots of corn 
crops; grassland. 

In a second stage, topographie features are considered from a digital élévation model 30 mètres 
grid size derivated from a 1/10000 scaled topographie map. Various topographie indexes for each 



sampling point are computed: the upslope cumulative area, the local slope and the Kirkby index 

(Kirkby, 1978) that détermine an occurence of soil saturation and defined as the logarithm of the ratio 

between the cumulative upslope area and the tangent of the local slope. Another topographie index 

referring to the downslope topography is also computed. It is an élévation différence, noted as DZ, 

and computed as follows: 
DZ jk = Zj-Zj k 

where Zj is the élévation at the data point and Zj k the élévation in the stream according to the flow 

pathway between the data point and the stream. In a first approximation, the lowest élévation of each 
axis is considered as Z\^. 

The water content data are averaged according to DZ classes of one mètre class interval to 
describe coarsely space and time distribution and model it. Around ten values are grouped in each 
class excepted for the last survey for which only half are considered. 

RESULTS 

Bulk data 

The spatial distribution of soil surface water content was analysed for each of the four surveys. The 
results on the bulk data are only presented for one survey and one depth, 27 January and 5-10 cm, 
because the same trends are observed for ail the surveys. 

Distinction of végétation cover brings no statification of the data, except for grassland which are 
mainly located in the valley bottom with high water content (Fig. 2A). One ont of the three soil types 
is well related to the water content data: the most wateriogged soils, mainly located along the 
channels, correspond to the greatest and the most scattered water content values (Fig. 2B). 

The water content data are not well related to the Kirkby's index (Fig. 3A): high water content, 
i.e. saturated conditions, does not always correspond to high index values. The local and upslope 
topography cannot explain the spatial distribution of the surface water content for ail data points. The 
data of water content are better related to the élévation différence to the stream DZ (Fig.3B). Two 
groups of values are distinguished: a group centered around an average value of 32% corresponds to 
DZ values higher than 4; a group distributed on the whole range of moisture corresponds to DZ 
values lower than 4. 

Data averaging and modelling 

The water content data averaged according to DZ classes of one mètre class interval confirms 
that there is a threshold beyond which water content is fairly constant and présents a low variability, 
and below which water content increases strongly and présents a higher variability as DZ decreases 
(Fig. 4). This threshold is mainly defined by the range of standard déviation that varies abruptly while 
the variations of mean values are both continuous and chaotic. This threshold corresponds to a mean 
topographie limit delineating the hillslope domain from the bottom domain. 



An exponential function defined as follows: 
W = A exp (-DZ / 2) + B 

is fitted (Fig. 4). It présents différent advantages: the exponential function is an often used synthetic 
draining function; it expresses a great then progressive^ slighter and siighter decrease of water 
content as DZ increases; it takes into account the iink between the values of the two previously 
defined domains. At the opposite it does not take into account the existence of a threshold that it is 
better defined by variations on standard déviations than on means. 

The map of the model residuals (Fig. 5), i.e. différence between measurements and model 
estimâtes for each data point, shows that the error sign is not randomly distributed in the catchment: 
in some areas, the errors are systematically positive or négative, located above ail in the bottom part 
and partly related to the location of the secondary drainage axis. This assumption is tested by taking 
into account a more précise drainage network derivated from a 1/1000 scaled topographie analysis 
used to calculate new values of DZ (Fig. 6). The variations of the water content are smoother, the 
slope at the origine is higher and the threshold is lower as a detailed topographie analysis is used 
(Fig. 7): the limit between the hillslope domain and the bottom domain is better defined. This 
calculation of DZ will be used in the further calculations. Conversely, if the élévation at the outlet of 
the catchment is chosen to calculate DZ (Fig. 6), according to a lower topographie accuracy or above 
ail to better take into account the emptying of the catchment in the driest conditions, the variations of 
the water content are more chaotic, the slope at the origine lower and the threshold higher (Fig. 7): 
the limit between the hillslope domain and the bottom domain of the catchment is more graduai and 
lesser defined even for the driest survey. 

The comparaison ofthe foursurveys 

The results observed for each survey and each depth are similar (Fig. 8 et 9). The élévation 
différence is relevant to distinguish the two domains on the catchment whatever the moisture status 
is: one of fairly constant water content that characterizes the hillslope domain, and the other of high 
and variable water content that characterize the bottom domain. The water content is approximately 
the same on hillslope for the two wettest surveys while it is lower for the two others: a mean soil 
surface water content on hillslope defines well the moisture status of the catchment. The curves are 
not parallel so that the Iink between the moisture of the two domains is loose. 

The study of variations on standard déviations does not show clearly if the threshold moves as 
the moisture on hillslope varies. It seems to be the same for the three wettest surveys and so must 
be considered as a fairly permanent delineation. It is uncertain for the driest survey for which the 
standard déviations are as high for hillslope as for the bottom. This is due, first, to the number of 
observations being halved; second, to the différence of évapotranspiration induced by the différent 
végétation covers being greater during summerthan during other periods. 

The decreasing exponential function was fitted to each of the four surveys, the two depths and 
the three tested élévation différence (Fig. 8 and 9). The fit obtained was always satisfactory (Fig. 8). 



The A and B parameters and the estimated accuracy (Tab. 1) were computed. B, that corresponds to 
a mean water content on the hiilslope domain, is more accurately determined than A. B is a good 
indicator of the variations in time of the moisture status on the catchment. A is more difficult to 
interpretate because it intégrâtes quick variations in space and time of moisture existing along the 
stream. The consistency of the results is best at 5-10 cm in depth and when taking into account a 
précise drainage network. 

DISCUSSION 

Data analysis 

AH the data analysis was canried out on gravimétrie water content measurements, although 
volumetric data were more relevant to analyse the variations of water storage on the catchment. This 
method was preferred because it allows to sample ail the catchment very quickly by five groups of 
people and the measurements really corresponds to the same moisture status. Measurements of 
volumetric water content by a cylindric method would take too long time while we had not a sufficient 
number of TDR that are moreover difficult in using in stone soils. To verify the consistency of the 
relationships established with the gravimétrie data, the bulk density measurements were used to 
calculate volumetric water contents at depth of 5-1 0 cm: the decrease of the water content as the 
élévation différence increases is higher and the two domains of the catchment are better defined 
when using volumetric than with gravimétrie data. In this context, this analysis is available even if 
using gravimétrie data. 

Physical significance 

Such spatial distributions of the soil surface water content, particularly the fairly constant 
moisture on the hiilslope domain have been already observed in this context by using a geostatistical 
approach (Gascuel-Odoux, 1984). The influence of the élévation différence to the stream on the 
spatial distribution of soil surface water content in the bottom domain is new. This élévation 
différence is an hydraulic head at least in saturated conditions. It may be considered as an effective 
property that can globally characterize the catchment but not point variations that are too erratic as it 
is shown by the analysis of the residuals from the model. The influence of this élévation différence 
may be interpreted not only according to topographie considérations but also to the spatial 
distribution of soils. 

As a matter of fact, as using the DEM for hydrological processes the Kirkby's index that referred 
to the upslope and the local topography is often used to explain space and time variations of soil 
moisture. Two arguments may explain that the élévation différence referring to downslope conditions 
is preferred here. Firstly the measurements are surface measurements and not vertical means of soil 
moisture: in the bottom domain where saturated conditions occur there is no différence; in the upper 
domain where the level of the water table differ in time and space the soil surface water content is 
fairly constant because the soils are the well drained (Gascuel-Odoux and Mérot, 1986), as the 



contrary to the vertical mean of soil moisture that présents continuous space and time variations. 
Secondly the measurements are realized out of the flood periods for which the downslope topography 
occurs more than the upslope topography in the draining proccesses. 

The soil distribution must be also considered because of the constant threshold on élévation 
différence regardless to the moisture status of the catchment and of the good relationship between 
surface moisture content and soil type: the bottom domain of the catchment may be seen as a 
permanent delineation and so partly due to intrinsic factors such as soil physical properties. This is 
apparently antagonic to the concept of variable contributive source area for which the surface 
saturation is highly time and space dépendant according to rainfall events. Indeed the extension of 
saturated areas varies with rainfall events on this site (Merot et al., 1993), but despite thèse 
variations the moisture remains relatively stable. It may be partly explained the high contrast of 
saturated hydraulic conductivity an order of magnitude higher for well drained slope soils as 
compared to the waterlogged bottom soils. 

Obviously expérimental validations must be realized to test the real interest of this index for soil 
distribution or soil surface moisture prédiction according to the geomorphological context and the 
accuracy of the DEM. 

Practical implications ofthe results 

The distinction of two levels of variability in the catchment has practical implications on sampling 
stratégies to evaluate the distribution of soil surface moisture on the catchment. Lower sampling 
densities must be chosen in the upper domain than in the bottom domain of catchments. This is 
important since the upper domain is much larger than the bottom domain. Moreover, field 
investigations limited to the upper domain of the catchment may be sufficient to assess the moisture 
status of the catchment. Thèse results also provide some insight on how to calibrate the information 
coming from a spacebome synthetic aperture imaging radar (Mérot et al., 1993). The data treatment 
on imaging radar will improve in using Windows to distinguish the hillslope from the bottom to 
compute the mean and the standard déviation on the radar signal. In a first phase, the élévation 
différence is a good index for this distinction as the DEM is available. 
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(Figures and tables) 



Figure 1 - Sampling scheme and location of the studied catchment within the Coet-Dan catchment 
(Brittany, France). 

Figure 2 - Water content on 27 January 1993 at 5-10 cm versus élévation différence according to: 
(A) végétation cover; (B) level of soil wateriogging. 

Figure 3 - Water content on 27 January 1993 at 5-10 cm versus: (A) Kirkby index; (B) élévation 
différence to the stream. 

Figure 4 - Mean (<W>), the standard déviation (<W>+sig, <W>-sig) and the fitted decreasing 
exponential function (<W> fitted) for water content on 27 January 1993 at 5-10 cm averaged by one 
mètre class intervais versus élévation différence to the stream. 

Figure 5 - Map of the residuals, expérimental data minus fitted values for each sample point, on 27 
January 1993 at 5-10 cm : (O) positive values; $S) négative values. 

Figure 6 - Schéma of the three methods using to calculate the élévation différence DZ from the 
stream. A: Z )A is located at the end of each axis, i.e. the main stream of the catchment; B: Z, B is 
located according to a drainage network derivated from a 1/1000 scaled topographie analysis; C: Z )C 
is the outlet of the catchment. 

Figure 7 - Mean water content on 27 January 1993 at 5-10 cm, averaged one mètre class intervais of 
élévation différence, versus élévation différence computed in référence to: (a) the outlet of the 
catchment; (b) the main stream of the catchment; (c) the secondary talwegs of the catchment. 

Figure 8 - Mean and standard déviation of water content of the four dates, at 0-5 cm, averaged one 
mètre class intervais, versus élévation différence computed from the secondary talwegs of the 
catchment. 

Figure 9 - Mean water content of the four dates, averaged one mètre class intervais, versus élévation 
différence computed from the main stream of the catchment at (A): 0-5 cm and (B): 5-10 cm. 

Table 1 - A and B parameters of the fitted exponential function: estimâtes and estimated standard 
déviation. 
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parameter A 


mean 


standard déviation 


A référence 


B référence 


C référence 


A référence 


B référence 


C référence 




02/18/1992 


18.6 


23.9 


87.6 


2.1 


1.2 


13.1 


0 - 5 cm 


10/24/1993 


34.9 


32.2 


49.6 


1.9 


5.3 


10.3 




01/27/1993 


36.8 


28.9 


53.4 


2.2 


8.1 


14.6 




09/03/1993 


10.8 


11.2 


12.3 


2.1 


1.8 


5.2 



5 - 10 cm 


10/24/1993 


20.4 


19.9 


34.1 


1.4 


4 


7.1 




01/27/1993 


28.9 


24.8 


39.2 


2.1 


5.3 


11.3 




09/03/1993 


7.5 


8.2 


1.9 


2.2 


2.2 


4.3 



parameter B 


mean 


standard déviation 


A référence 


B référence 


C référence 


A référence 


B référence 


C référence 




02/18/1992 


31.7 


32.3 


37.2 


0.7 


0.4 


1.4 


0 - 5 cm 


10/24/1993 


23.2 


26.8 


29.3 


0.6 


1.8 


1.7 




01/27/1993 


32.5 


36.7 


40.5 


0.8 


2.7 


2.5 




09/03/1993 


7.9 


8.3 


9.9 


0.5 


0.5 


1 




5 - 10 cm 


10/24/1993 


25.2 


27.2 


28 


0.5 


1.4 


1.2 




01/27/1993 


30.7 


33.25 


37.3 


0.7 


1.8 


1.9 




09/03/1993 


14.7 


14.7 


16.7 


0.5 


0.5 


0.8 



table 1 



3. conclusion 

L'étude expérimentale décrite dans ce chapitre avait pour objectif premier, d'utiliser les 
images ERS1 pour étudier l'évolution de l'humidité de surface sur des bassins versants 
d'ordre 1 et en déduire la structure du réseau hydrographique élargi aux zones 
contributives. La complexité des images radar n'autorise pas une exploitation simple de 
l'information radar pour calculer l'humidité des sols. Les mesures effectuées dans le 
cadre d'une validation terrain des images radar ont cependant permis de distinguer deux 
systèmes: un versant et un bas-fond. Le premier est caractérisé par une humidité 
pondérale homogène. Le second présente en revanche des valeurs d'humidité très 
variables et supérieures à la capacité au champ. L'interface entre ces deux systèmes est 
délimitée par un seuil de potentiel référencé par rapport aux axes de drainages principaux. 

Malheureusement il n'a pas été possible d'établir avec certitude la stabilité ou 
l'évolution du seuil en fonction du degré global d'humidité du bassin versant, et de 
proposer une interprétation physique en terme de couplage entre le réseau et les zones 
saturées. D'autres études expérimentales seront donc nécessaires pour définir la 
dynamique du seuil en fonction du référentiel aval pendant un cycle saisonnier complet. 



189 



CONCLUSION 
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Ce travail s'inscrit dans le cadre d'un apport à la compréhension de l'organisation des 
réseaux hydrographiques et du relief. Deux approches complémentaires ont permis 
d'établir certains traits fondamentaux des ces deux systèmes évolutifs. 

La première approche est une analyse statistique des caractéristiques géométriques de 
deux fleuves situés en Bretagne : le Blavet et la Vilaine. Elle montre l'existence de deux 
types de distributions de longueurs de cours d'eau et d'aire drainée. On observe: (i) des 
lois exponentielles des longueurs de cours d'eau entre les confluences, qui traduisent une 
distribution aléatoire des confluences au sein de l'arborescence hydrographique, (ii) des 
distributions de longueurs cumulées et d'aire drainée qui suivent des lois de puissances 
caractéristiques des systèmes d'agrégation aléatoire gouvernés par une contrainte de 
branchement vers l'aval. Les propriétés géométriques des réseaux hydrographiques 
vérifient les lois des arborescences au sens large, telles que les lois de Horton. Mais ces 
dernières ne représentent en rien une spécificité des système hydrographiques et ne 
peuvent de ce fait être utilisée comme critère de validation dans les modèles de 
configuration de réseaux hydrographiques. Nous observons également une remarquable 
homogénéité de la densité de cours d'eau en Bretagne. Ce résultat est en contradiction 
avec les récents modèles d'organisation fractale. Il peut suggérer une modification des 
propriétés géométriques des réseaux au cours de leur croissance. Une compréhension 
complète de l'organisation spatiale des réseaux hydrographiques passe donc par l'étude 
de leur dynamique. 

Les distributions statistiques de longueurs (cumulées ou non) des deux fleuves bretons 
étudiés, ne suivent les lois de régressions décrites plus haut qu'au-delà d'une certaine 
longueur seuil de l'ordre de 1 km qui correspond à la longueur moyenne des biefs 
d'ordre 1. Ce seuil est à la fois une limite inférieure à l'extension des réseaux 
hydrographiques et une longueur caractéristique du relief. 

La deuxième approche vise à établir un modèle de croissance des réseaux 
hydrographiques afin de déterminer l'influence des processus physiques fondamentaux 
sur l'évolution conjuguée des réseaux et du relief. Le modèle proposé est construit à la 
fois sur la base de règles déterministes et stochastiques. Le déplacement de volume d'eau 
est simulé par le déplacement de marcheurs élémentaires (précipitons) qui évoluent sur la 
topographie tout en la modifiant par des règles simples d'érosion/sédimentation. On 
traduit ainsi de manière simple, en terme de règles de conduites de marcheurs, les 
principaux processus physiques d'érosion et de sédimentation (diffusion, transport 
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advectif, ...)• Comme la plupart des modèles récents, nous pensons que la dynamique 
des systèmes géomorphologiques est fortement contrôlé par une modification des règles 
d'érosion lorsqu'on passe du versant (écoulement diffus souvent en milieu poreux) au 
réseau hydrographiqe (écoulement à surface libre). L'originalité de notre simulation est 

(i) de ne pas exiger a priori de conditions de branchement aval des réseaux, et (ii) 
d'utiliser une condition sur la valeur du flux hydraulique pour déterminer le régime des 
écoulements hydrauliques. 

Avec de telles hypothèses, nous montrons plusieurs modes de croissance des réseaux 
hydrographiques : (i) une croissance vers l'amont à partir des limites les plus basses du 
système, (ii) une concentration des flux transversalement à la pente moyenne, (iii) des 
phénomènes de capture de sous-réseaux établis. La prédominance d'un mode par rapport 
à l'autre dépend beaucoup de l'amplitude de la perturbation imposée aux limites du 
système. Si la perturbation est forte (abaissement du niveau de la mer, par exmple), le 
mode (i) sera dominant. Sinon, le réseau hydrographique s'organisera suivant les modes 

(ii) et (iii). 

Nous avons aussi étudié de façon quantitative le développement d'instabilités de pente 
lié aux différents mécanismes d'érosion. Les simulations retrouvent les résultats proposés 
par les études d'instabilités marginales de Smith et Bretherton (1972) ou Loewenhertz 
(1991). Les instabilités ne peuvent se développer que s'il existe une relation non-linéaire 
entre le flux de matière et le flux hydraulique. C'est le cas des simulations qui prennent en 
compte une modification du régime hydraulique de type versant-rivière. Les instabilités se 
développent alors perpendiculairement à la pente principale avec une facteur 
d'amplification qui dépend des paramètres de la loi d'érosion. 

Dans la troisième partie de ce travail, par une approche expérimentale, nous avons 
chercher à établir des critères qui permettraient de localiser cette transition versant/réseau 
dans le cas des régions humides. Le rôle particulier des zones saturées dans la cinétique 
des transferts de volumes d'eau, nous a amené à utiliser leur extension géographique 
comme marqueur de la localisation de cette transition. 
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Effect of saturated areas on backscattering coefficient 
of the ERS 1 synthetic aperture radar: First results 
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Abstract. One way to estimate spatial distribution of water content at the soil surface 
consists of using active microwave remote sensing. It has been theoretically and 
experimentally demonstrated for unsaturated conditions. Neverthéless, radar data are 
ambiguous when ponding conditions occur, as in variable source areas, due to the 
contradictory influence of the dielectric effect and the specular effect on the 
backscattering atténuation coefficient, <x 0 . A procédure is considered, based on a 
topographie analysis, to take into account the influence of the two effects on radar 
response. The présent results stress the capabilities of the ERS (European Remote 
Sensing Satellite) radar to survey saturated areas in time and space. 



Introduction 

Mapping areas where an excess of water occurs at the soil 
surface is of major importance in hydrology. Surface runoff 
in the widest sensé mày be generated by two main mecha- 
nisms. First, Hortonian overland flow ("infiltration excess 
runoff") occurs when the rainfall intensity exceeds the soil 
infiltration capacity. This is independent of the watershed 
morphology and of position in the landscape; for homoge- 
neous soils it involves the whole watershed considered. 
Second, saturation overland flow and retum flow, according 
to the variable source area concept, which was first pro- 
posed by Cappus [1960] and then developed by différent 
authors [Kirkby and Chorley, 1967; Beven and Wood, 1983; 
Merot, 1988] represent surface runoff generated by rainfall 
on saturated areas. According to this concept, not ail the 
watershed area will contribute equally to runoff; the runoff- 
generating zone may be of varying size and represents the 
ground water exfiltration areas. Its occurrence dépends both 
on topographical factors, specifically convergence of runoff 
flow and slope and on soil characteristics, specifically per- 
meability. The distribution and the extent of thèse areas 
control the saturation overland flow and return flow. Différ- 
ent authors [Beven and Kirkby, 1979; Grayson et al., 1992] 
have developed models based on this second concept. 

The difficulty of mapping the runoff exfiltration zone, due 
to the rapid variation in space and time of thèse areas, has 
often been pointed out. It is the main limit to validating the 
variable source area concept and the corresponding models. 

One way to estimate water content at the soil surface 
consists of using active microwave remote sensing. This 
technique has been theoretically and experimentally demon- 
strated [Chanzy, 1991; Engman et al., 1989; Mo et al., 1984; 
Ulaby et al., 1978; Van de Griend and Engman, 1985]. 
However, when ponding conditions occur, as in variable 
source areas, radar data are ambiguous due to the contra- 
Copyright 1994 by the American Geophysical Union. 
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dictory influence of the dielectric effect and the specular 
effect on the backscattering atténuation coefficient, a 0 . The 
dielectric effect is caused by the dielectric capacity of soils, 
which increases with the water content; the backscattering 
coefficient has a positive corrélation with the volumetric 
water content of soils until a threshold close to field capac- 
ity. The specular effect is due to the reflection at the 
saturated soil surface or on ponding water that involves a 
decreasing backscattering coefficient when the saturated soil 
surface increases. 

In a previous paper, Brun et al. [1990] showed that 
saturated areas can be detected using radar data. Thèse data 
came from a helicopter-borne instrument, ERASME (Etude 
par Radar des Sols de la Mer) with the same configuration as 
the ERS 1 (European Remote Sensing) satellite synthetic 
aperture radar (SAR); under a threshold backscattering 
value the measured plots corresponded to saturated areas; 
over this threshold value they corresponded to soils at field 
capacity. In another paper [Merot and Chanzy, 1991] an 
experiment characterizing the SAR response in controlled 
ponding conditions was reported. The extent of surface 
saturation was measured on two plots with différent rough- 
nesses. The extent varied in time from 50% to 0%. The 
surface saturation was simultaneously measured by the 
crane-mounted radar RAMSES with two configurations (5.3 
GHz, VV, 23°, and 5.3 GHz, HH, 15°) and by photography. 
It was shown, on pixels of 10 m 2 size, that the backscattering 
coefficient decreases when the percentage of saturated area 
increases. The corrélation gradient was the same for the 
différent roughnesses. Because of ambiguous SAR values 
between measurements in ponding conditions and measure- 
ments in drier conditions, the necessity of complementary 
data was emphasized. 

In this paper, as a part of a European Space Agency (ESA) 
ERS 1 project [Louhala et al., 1993], preliminary results 
from ERS 1 SAR are presented. The objective is to analyze 
the radar data and to emphasize an interprétation procédure 
taking into account the two processes of backscattering 
atténuation when the soil is drying and ponding. 
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Materials and Methods 

The expérimental watershed is situated in the center of 
Brittany, in France, at the intersection of the descending and 
ascending swathes of ERS 1 during the "ice phase" period 
(Lambert coordinates: 215.8 and 2,343.79). The study water- 
shed was described by Brun et al. [1990] and C. Loumagne 
et al. (Evaluation of the ERS-1/SAR capacity to estimate 
surface soil moisture, submitted to Remote Sensing of 
Environment, 1993). It drains an area of 2 km 2 into the upper 
part of the Coët-Dan River, which is a tributary of the Evel 
River. The relief of the catchment area is smooth. The 
substratum consists of brioverian schists. Thèse silty schists 
alternate locally with sandstone. Alluvial material covers the 
bottom of the valley. Soils consist of loamy altération of 
schists, mixed with eolian loam on west hillslopes. The top 
of the slopes are overlaid by dystric or aquic eutrochrepts 
(brown acidic, weakly leached soils). Glossaqualfs (degrad- 
ed hydrorhorphic soils) and Auvents (alluvial soils) cover low 
grounds. The mean annual rainfall is 711 mm. Rainfall of 
over 20 mm per day is exceptional. The basin is character- 
ized by intensive farming, with maize, wheat, temporary and 
permanent pasture for milk production, indoor stock farm- 
ing, and vegetables. 

Saturated areas spread along the river. Ground-based 
gravimétrie water content measurements were carried out 
during différent campaigns simultaneously with the SAR 
satellite measurements. Field measurements were made 
along 10 transects from the top to the bottom of hillslopes, 
with a 25 m spacing, at 0-5 and 5-10 cm depths (about 500 
samples). An automatic rain gage located at the top of the 
watershed gave the rainfall amount. Land use was surveyed 
during the campaign. It mainly includes maize (after harvest- 
ing), winter wheat, and pasture. 

Radar measurements were made every 3 day s, except for 



two scènes that are missing, from January to March 1992 
(the so-called "ice phase") and every 35 days from April to 
December 1992. A preliminary processing of the scènes was 
made by ESA at the ESRIN (European Space Research 
Institute, Frascati, Italy). The pixel size is 12.5 x 12.5 m . 

For the purpose of the study, contiguous cells of approx- 
imately the same size (about 1 or 2 ha), with various soils and 
végétation, were selected from the top to the bottom of each 
hillslope where ground-based water content was measured 
(Figure 1). Soil and végétation types of the cells are reported 
in Table 1. The mean of the backscattering values of the 
différent pixels (between 40 and 70 pixels, depending on the 
size of the cell) of a cell was attributed to the cell. 

The short time lag between two images during the "ice 
phase" makes it possible to compare the value of the 
backscattering values from date to date; the only parameter 
susceptible to rapid variation is soil water content. Therefore 
the analysis will be based on the variation of the backscat- 
tering values between one day and another, and not on the 
absolute value. 

During the ice phase, the weather was unfortunately not 
favorable to the experiment; the beginning of the period was 
cold, and soil surface frost occurred during satellite mea- 
surements, on day 28 and day 52. Moreover, only two rainy 
periods occurred, from day 37 to 46 and, at the end, from day 
70 to 80 (Figure 2). 



Results 

In Figure 3 the absolute backscattering values are plotted 
for différent cells. There is a large range of values for the 
différent cells at each date. Two groups of cells are distin- 
guished: The most numerous one (cells 61, 63, 65, 67, 74, 75, 
76. 77) is composed of cells having simultaneous variations 
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Table 1. 


Soil Type and Végétation Cover During Winter 1992 for Ail Cells 






Type of 


n UlilUC-1 




Végétation 


of Cells 


Type of Soil 


Cover 


61, 62, 63 


dystric or aquic eutrochrept, more than 80 cm 


winter wheat 


67 


dystric or aquic eutrochrept, more than 80 cm 


corn 


74 


dystric or aquic eutrochrept, 60-80 cm 


meadow 


75, 76, 77 


dystric or aquic eutrochrept, 60-80 cm 


plowed field 


69 


glossaqualf on clayey alterites 


winter wheat 


70, 71 


giossaqualf on clayey alterites 


vegetable 


72 


glossaqualf on clayey alterites 


meadow 



in time of the backscattering values. The other group (cells 
69, 70, 71) has a more heterogeneous radar response and 
lower values and présents smoother variations in time, or 
variations that are not cohérent with the previous one. This 
allows us to divide the population of cells in two. 

In Figure 4 the différences between the backscattering 
values for two consécutive observations, from day d to day 
d + 3 (when scènes are missing day d + 6 is used), are 
plotted for différent cells; each of the two groups of cells is 
presented separately. Cell 72 was added to the second group 
in spite of its intermediary behavior. Comparing to Figure 2, 
the first group (Figure 4a) includes cells with backscattering 
values that increase when rainfall occurs (i.e., day 40), and 
decrease afterward, in contrast to the second one (Figure 
4b). Backscattering values also strongly decrease when frost 
occurs (day 52), due tô the low dielectric constant of the ice. 

The two groups are not randomly distributed in the 
landscape. None of the cells of the first group is situated in 
the bottom of the talweg. On the contrary, cells of the 
second group are at the foot of the hillslopes or in the bottom 
of the talweg. Field observations show particularly strong 
variations in time of saturated areas in cell 71, which is 
located in the head of the saturated talweg. 

In Figure 5, the behavior of one cell from each of the two 
groups is plotted; cell 76 is at the top of a transect, and cell 
71 is at the bottom. The two cells have an inverse behavior, 
showing that the backscattering coefficient increases for one 
cell, whereas it decreases for the other. 

Interprétation 

The main features of the spatial distribution of soil water 
content in such a watershed during winter were briefly 
described in the introduction. Saturated areas are located 
near the river (the variable source area), and homogeneous 
values of water content close to field capacity are observed 
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on the hillslope. Figure 6 présents a comprehensive diagram 
of the gravimétrie water content in the soil surface plotted 
against the relative élévation from the river. Again, two 
groups occur: At the top and on the hillslope, water content 
is moderate and of small standard déviation, whereas the 
second group, whose relative élévation from the talweg axis 
is less than 5 m, shows a significant gradient and a large 
standard déviation, partly corresponding to saturated areas. 
The large standard déviation is due to the heterogeneous 
distribution of ponding water, depending on soil roughness. 

The following scheme is proposed to explain the twofold 
behavior of the SAR signal. When the target of the SAR is 
the bottom of the talweg where the soils are always saturated 
in winter, the backscattering value is low (due to the 
specular effect) and has few variations (Figure 4b, cells 69 
and 70). When the target of the SAR is just upward from the 
talweg in the zone of source area variations, the backscat- 
tering value decreases when it rains (again due to the 
specular effect), and increases afterward, due to the varia- 
tion of the total surface of the saturated area [Merot and 
Chanzy, 1991]. When the target of the SAR is on the slope in 
an upper zone where no water excess occurs, the backscat- 
tering value increases when it rains (due to the dielectric 
effect). The variations of the backscattering values along a 
transect from the talweg to the top (cells from 70 to 77) 
during a short rainy period (between day 37 and day 46) are 
presented in Figure 7. Variations along the hillslope corre- 
spond to the scheme described here. 

Another argument that confirms the effect of the saturated 
zone cornes from analysis of the same two plots 71 and 76 
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Figure 2. Rainfall. 
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Figure 3. Absolute backscattering values for différent cells 
versus time. 
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Label cells 





61 






63 




65 






67 




74 






75 






76 






77 




mean and standard déviation 



Figure 4. Différences between day /' and day i - 3 in mean 
backscattering values versus time (a) for cells corresponding 
to two hillslopes (south and north); (b) for cells correspond- 
ing to the talweg. 



during the whole experiment period, including sutnmer (Fig- 
ure 8): The two plots show an opposite variation in winter, 
when saturated areas occur on plot 71; then, in summer, 
when saturated areas disappear, the curves cross over and 
show cohérent variations with the dielectric effect. 

A last comment has to be made on using radar data with a 
temporary soil frost: As it was shown, the soil surface was 
frozen from time to time, because of the unexpected cold 
weather. Consequently, backscattering values strongly de- 
creased, impeding any interprétation in terms of soil mois- 
ture or saturated areas with such weather conditions. 
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Figure 5. Différences in mean backscattering value versus 
time for cell 76 (top of a transect) and 71 (bottom). 
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relative élévation from river in m 



Figure 6. Comprehensive diagram of the gravimétrie water 
content in surface soil (0-5 cm depth) versus the relative 
élévation from the river on two dates. 



Conclusion 

Thèse preliminary results stress the capabilities of ERS 1 
radar to survey saturated areas. They emphasize the impact 
of saturated areas on the backscattering coefficient, gener- 
ally underestimated or neglected. A possible explanation of 
the opposite variations of the backscattering value observed 
after a rainfall, due to two contradictory physical effects, is 
proposed. We have to take into account the position in the 
landscape of the target cells, and to have in mind a concep- 
tual scheme of the spatial distribution of saturated areas in 
the landscape. A topographie modeling has to be done, to 
sort the SAR data into two groups, the fîrst one depending on 
the dielectric effect and the second on the specular effect. 
From a hydrological point of view, ERS 1 radar seems to be 
a relevant tool for validation of the variable source area 
concept and the corresponding models. 




73 74 75 76 

Number of cells on transect 70 - 77 



Figure 7. Différences in mean backscattering value during 
a rainfall period from the top to the bottom of a hillslope: cr 0 
between day 37 and days 40, 43, 46. 
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Figure 8. Mean backscattering values versus time for cells 
76 and 77 (top of a transect) and cell 71 (bottom) for the 
whole measurement period. 
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Les réseaux hydrographiques jouent un rôle prépondérant dans la régulation des 
flux d'eau et de matière à la surface des continents. Sièges d'une érosion intense, ils 
contrebalancent, en découpant le paysage par un réseau de vallées, les apports 
endogènes d'origine tectonique. Leur structure et leur degré d'extension et de 
ramification sont des paramètres clés dans l'évolution du relief sur des pas de temps 
longs (10- -10 7 années), et dans la propagation d'ondes de crues sur des pas de 
temps courts (1 heure, qques jours). Comprendre leur évolution est un enjeu 
intéressant à la fois d'un point de vue fondamental, pour la détermination des 
processus physiques fondamentaux de transferts de matière à des échelles variant 
du continent au versant, mais également dans une perspective d'application à des 
besoins de prévision hydrologiques. 

Le travail présenté s'attache d'une part à définir un schéma d'organisation 
possible et, d'autre part, à définir un modèle simple d'extension et d'évolution des 
réseaux hydrographiques et du relief. L'analyse statistique des propriétés 
géométriques de deux réseaux hydrographiques bretons : la Vilaine et le Blavet, met 
en évidence, entre autre, une distribution aléatoire des confluences au sein du 
système hydrographique et une densité homogène sur 4 ordres de grandeur 
d'échelle. Ces deux observations témoignent d'une géométrie aléatoire et oriente la 
modélisation vers un modèle mixte stochastique-déterministe. Le 
modèle proposé traduit la physique des principaux processus ========== 

d'érosion par des règles simples de déplacement et d'action de " 
marcheurs lancés aléatoirement sur une topographie. D'utilisation 
très souple, cet outil offre la possibilité d'aborder l'étude de ~^ 
l'évolution des reliefs et des axes de drainage en fonction de -- 
l'importance relative des processus fondamentaux que sont les 
déplacements et les transports de matière par diffusion et advection. - ■ - 
il est ainsi possible de simuler la croissance de perturbations . 
suivant la prépondérance de processus particulier. 
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Réseaux hydrographiques, lois d'échelle, modélisation, érosion. 
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